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The so called /(-R)-gravity has recently attracted a lot of interest since it could be, in principle, 
able to explain the accelerated expansion of the Universe without adding unknown forms of dark 
energy/dark matter but, more simply, extending the General Relativity by generic functions of the 
Ricci scalar. However, apart several phenomenological models, there is no final /(i?)-theory capable 
of fitting all the observations and addressing all the issues related to the presence of dark energy 
and dark matter. An alternative approach could be to "reconstruct" the form of f{R) starting from 
data without imposing particular classes of model. In this review paper, we will consider two typical 
cosmological problems where the role of dark energy and dark matter is crucial. Firstly, assuming 
generic /(i?), we show that it is possible to relate the cosmographic parameters (namely the deceler- 
ation go, the jerk jo, the snap so and the lerk lo parameters) to the present day values of f{R) and its 
derivatives /<"'(-R) = d"f/dR" (with n = 1,2,3) thus offering a new tool to constrain such higher 
order models. Our analysis gives the possibility to relate the model-independent results coming from 
cosmography to some theoretically motivated assumptions of f{R) cosmology. Besides, adopting 
the same philosophy, we take into account the possibility that galaxy cluster masses, estimated at 
X-ray wavelengths, could be explained, without dark matter, reconstructing the weak-field limit of 
analytic ,f{R) models. The corrected gravitational potential, obtained in this approximation, is used 
to estimate the total mass of a sample of 12 well-shaped clusters of galaxies. Results show that such 
a gravitational potential provides a fair fit to the mass of visible matter (i.e. gas + stars) estimated 
by X-ray observations, without the need of additional dark matter while the size of the clusters, as 
already observed at different scale for galaxies, strictly depends on the interaction lengths of the 
corrections to the Newtonian potential. These two examples could be paradigmatic to overcome 
dark energy and dark matter problems by the extended gravity approach. 



I. INTRODUCTION 



As soon as astrophysicists realized that Type la Supernovae (SNela) were standard candles, it appeared evident 
that their high luminosity should make it possible to build a Hubble diagram, i.e. a plot of the distance - redshift 
relation, over cosmologically interesting distance ranges. Motivated by this attractive consideration, two independent 
teams started SNela surveys leading to the unexpected discovery that th e Universe expansio n is speeding up rather 
than decelerating as assumed by the Cosmological Standard Model [zl, Il25 
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has now been strengthened by more recent data coming from SNela surveys [a]Ta.l60rr96l.ll33l. fisi . [iH^JlTi, large 
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l85l . Il21l . [12a, Il54l | and cosmic microwave background (CMBR) anisotropy spectrum 
. This large data set coherently points toward the picture of a spatially fiat Universe undergoing 
an accelerated expansion driven by a dominant negative pressure fluid, typically referred to as dark energy [50| . 

While there is a wide consensus on the above scenario depicted by such good quality data, there is a similarly wide 
range of contrasting propos als to solve the dark energy puzzle. Surprisingly, the sim ples t explanat ion, namely the 
cosmological constant A Sim, is also the best one from a statistical point of view |l45l . Il56l . Il57l | . Unfortunately, 
the well known coincidence and 120 orders of magnitude problems render A a rather unattractive s olution fr om a 
theoretical point of view. Inspired by the analogy with inflation, a scalar fleld 0, dubbed quintessence [ll9l . Il2d |. has 
then been proposed to give a dynamical A term in order to both flt the data and avoid the above problems. However, 
such models are still plagued by difficulties on their own, such as the almost complete freedom in the choice of the 
scalar field potential and the fine tuning of the initial conditions. Needless to say, a plethora of alternative models 
are now on the market all sharing the main property to be in agreement with observations, but relying on completely 
different physics. 

Notwithstanding their differences, all dark energy models assume that the observed apparent acceleration is the 
outcome of some unknown ingredient, at fundamental level, to be added to the cosmic pie. In terms of the Einstein 
equations, G^^ = x^^ii/, the right hand side should include something more than the usual matter and radiation 
components in the stress - energy tensor. 

As a radically different approach, one can also try to leave unchanged the source side (actually "observed" since 
composed by radiation and baryonic matter), but rather modifying the left hand side. In a sense, one is therefore 
interpreting cosmic speed up as a first signal of the breakdown of the laws of physics as described by the standard 
General Relativity (GR). Since this theory has been experimentally tested only up to the Solar System scale, there is 
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no a priori theoretical motivation to extend its validity to extraordinarily larger scales such as the extragalactic and 
cosmological ones (up to the last scattering surface!). Extending GR, not giving up to its positive results at lo cal s c ales , 
opens the way to a large class of alternative theo ries o f gravity ranging from extra - dimensions [gI, [13, [gI, llQll . Il02l | 
to nonminimally cou pled scalar fields [36, 61, "ts", ^12^ . In particular, we are interested here in fourth order theories 
0, HH, m, [H, m, m, m, [H, OM, JM, mj] based on replacing the scalar curvature R in the Hilbert-Einstein 

action with a generic analytic function f{R) which should be reconstructed starting from data and physically motivated 
issues. Also referred to as /(i?) -gravity, some of these models have been shown to be able to both fit the cosmological 
data and evade the Solar System constraints in several physically interesting cases BHQia, 151, 160]. 

In this review paper, we will face two of the main problems directly related to the dark energy and dark matter 
issues: cosmography and clusters of galaxies. These are typical examples where the standard General Relativity and 
Newtonian potential schemes fail to describe dynamics since data present accelerated expansion and missing matter. 
Our goal is to address them by /(i?)-gravity. 

A. Cosmography: why? 

It is worth noting that both dark energy models and modified gravity theories seem to be in agreement with 
data. As a consequence, unless higher precision probes of the expansion rate and the growth of structure will be 
available, these two rival approaches could not be discriminated. This confusion about the theoretical background 
suggests that a more conservative approach to the problem of cosmic acceleration, relying on as less model dependent 
quantities as possible, is welcome. A possible solution could be to come back to the cosmography |l7Cll | rather than 
finding out solutions of the Friedmann equations and testing them. Being only related to the derivatives of the scale 
factor, the cosmographic parameters make it possible to fit the data on the distance - redshift relation without any a 
priori assumption on the underlying cosmological model: in this case, the only assumption is that the metric is the 
Robertson - Walker one (and hence not relying on the solution of cosmological equations) . Almost eighty years after 
Hubble discovery of the expansion of the Universe, we can now extend, in principle, cosmography well beyond the 
search for the value of the only Hubble constant. The SNela Hubble diagram extends up to z = 1.7 thus invoking the 
need for, at least, a fifth order Taylor expansion of the scale factor in order to give a reliable approximation of the 
distance - redshift relation. As a consequence, it could be, in principle, possible to estimate up to five cosmographic 
parameters, although the still too small data set available does not allow to get a precise and realistic determination 
of all of them. 

Once these quantities have been determined, one could use them to put constraints on the models. In a sense, 
we can revert the usual approach, consisting in deriving the cosmographic parameters as a sort of byproduct of an 
assumed theory. Here, we follow the other way around expressing the model characterizing quantities as a function of 
the cosmographic parameters. Such a program is particularly suited for the study of fourth order theories of gravity. 
As it is well known, the mathematical difficulties entering the solution of fourth order field equations make it quite 
problematic to find out analytical expressions for the scale factor and hence predict the values of the cosmographic 
parameters. A key role in /(i?)-gravity is played by the choice of the f{R) function. Under quite general hypotheses, 
we will derive useful relations among the cosmographic parameters and the present day value of f^^\R) = cPf/dR" , 
with n = 0, . . . , 3, whatever ,f{R) is^. Once the cosmographic parameters will be determined, this method will allow 
us to investigate the cosmography of f{R) theories. 

It is worth stressing that the definition of the cosmographic parameters only relies on the assumption of the 
Robertson - Walker metric. As such, it is however difficult to state a priori to what extent the fifth order expansion 
provides an accurate enough description of the quantities of interest. Actually, the number of cosmographic parameters 
to be used depends on the problem one is interested in. As we will see later, we are here concerned only with the SNela 
Hubble diagram so that we have to check that the distance modulus ficpiz) obtained using the fifth order expansion 
of the scale factor is the same (within the errors) as the one ^de{z) of the underlying physical model. Being such 
a model of course unknown, one can adopt a phenomenological parameterization for the dark energy^ equation of 
state (EoS) and look at the percentage deviation A/i//i£)£ as function of the EoS parameters. We have carried out 
such exercise using the CPL model, introduced below, and verified that IS.ii/ijlde is an increasing function of z (as 
expected), but still remains smaller than 2% up to z ~ 2 over a wide range of the CPL parameter space. On the other 



^ As an im portant remark, we stress that our derivation will rely on the metric formulation of f(R) theories, while we refer the reader to 

[l27l .ll28l| for a similar work in the Palatini approach. 
^ Note that one can always use a phenomenological dark energy model to get a reliable estimate of the scale factor evolution even if the 

correct model is a fourth order one. 
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hand, halting the Taylor expansion to a lower order may introduce significant deviation for z > 1 that can potentially 
bias the analysis if the measurement errors are as small as those predicted by future SNela surveys. We are therefore 
confident that our fifth order expansion is both sufficient to get an accurate distance modulus over the redshift range 
probed by SNela and necessary to avoid dangerous biases. 



In the second part of this review we will apply the /(i?)-gravity approach to cluster of galaxies. In fact, changing 
the gravity sector has consequences not only at cosmological scales, but also at galactic and cluster scales so that 
it is mandatory to investi gate the low energy limit of su ch th eories. A strong debate is open with different results 
arguing in favor [sl. [35l. I4ll. l63 . Iiod Il47l | or against [3, 0, such models at local scales. It is worth noting that, as 
a general resul t, higher order theo ries of gravity cause the gravitational potential to deviate from its Newtonian 1 jr 
scali ng lisl . losl . Il05l . ll42l . Il46l . Il52l | even if such deviations may be vanishing. 

In [2g|, the Newtonian limit of power law /(i?) — /o-R" theories has been investigated, assuming that the metric in 
the low energy limit (<I>/c^ << 1) may be taken as Schwarzschild- like. It turns out that a power law term {r/rc)^ 
has to be added to the Newtonian 1/r term in order to get the correct gravitational potential. While the parameter 
(3 may be expressed analytically as a function of the slope n of the f{R) theory, Tc sets the scale where the correction 
term starts being significant. A particular range of values of n has been investigated so that the corrective term is an 
increasing function of the radius r thus causing an increase of the rotation curve with respect to the Newtonian one 
and offering the possibility to fit the galaxy rotation curves without the need of further dark matter components. 

A set of low surface brightness (LSB) galaxies with extended and well measured rotation curves has been considered 
[H^-IsbI. These systems are supposed to be dark matter dominated, and successfully fitting data without dark matter is 
a strong evidence in favor of the approach (see also [tHI for an independent analysis using another sample of galaxies) . 
Combined with the hints coming from the cosmological applications, one should have, in principle, the possibility to 
address bot h th e dark energy and dark matter problems resorting to the same well motivated fundamental theory 
[iB, [13, llOd l. Nevertheless, the simple power law f{R) gravity is nothing else but a toy-model which fail if one 
tries to achieve a comprehensive model for all the cosmological dynamics, ranging from the early Universe, to the 
large scale structure up to the late accelerated era jH, [o^l . 

A fundamental issue is related to clusters and superclusters of galaxies. Such structures, essentially, rule the large 
scale structure, and are the intermediate step between galaxies and cosmology. As the galaxies, they appear dark- 
matter dominated but the distribution of dark matter component seems clustered and organized in a very different 
way with respect to galaxies. It seems that dark matter is ruled by the scale and also its fundamental nature could 
depend on the scale. For a comprehensive review see Q. 

In the philosophy of /(i?) -gravity, the issue is to reconstruct the mass profile of clusters without dark matter, i.e. 
to find out corrections to the Newton potential producing the same dynamics as dark matter but starting from a well 
motivated theory. 

In conclusion, /(i?)-gravity, as the simplest approach to any extended or alternative gravity scheme, could be the 
paradigm to interpret dark energy and dark matter as curvature effects acting at scales larger than those where 
General Relativity has been actually investigated and probed. 

Let us discuss now how cosmography and then galaxy clusters could be two main examples to realize this program. 



The key rule in cosmography is the Taylor series expansion of the scale factor with respect to the cosmic time. To 
this aim, it is convenient to introduce the following functions: 



B. Clusters of galaxies: why? 



II. THE COSMOGRAPHIC APPARATUS 
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which are usually referred to as the Hubble, deceleration, jerk, snap and lerk parameters, respectively. It is then a 
matter of algebra to demonstrate the following useful relations : 



H=H^{j + iq + 2) , 
dPH/dt^^H^ [s - 4j - iq{q + 4) - 6] , 
d^H/dt'^=H'' [l-5s + 10{q + 2)j + 30{q + 2)q + 24] 



(6) 
(7) 
(8) 
(9) 



where a dot denotes derivative with respect to the cosmic time t. Eqs.(l6]) - ((9| make it possible to relate the derivative 
of the Hubble parameter to the other cosmographic parameters. The distance - redshift relation may then be obtained 
starting from the Taylor expansion of a{t) along the lines described in [iol . Il66l . Il69| . 



A. The scale-factor series 



With these definitions the series expansion to the 5th order in time of the scale factor will be: 

a{t)=a{t,) Uo{t to) ^Hl{t - t^f + ^^H',{t - tof + '^H^it to)' + ^-^H^oit " ^o)' + 0[{t - tof\ \ (10) 
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^2 I Jo jjZ 



So 



lo 



a{to) — ""^^ 2 ' 3! ' 4! '5! 

It's easy to see that Eq. (fTT|) is the inverse of redshift z, being the redshift defined by: 
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The physical distance travelled by a photon that is emitted at time t* and absorbed at the current epoch to is 

D = c j dt^ c{to - i*) 

Assuming t^ — to — ^ and inserting in Eq. (fTT|) we have: 

a{to) 1 . X 

l + Z= "^n, = —7—. (12) 
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The inverse of this expression will be: 
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Then we reverse the series z{D) D{z) to have the physical distance D expressed as function of redshift z: 
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From this we have: 



D{z) = ^{Vl+V\z + Vlz' + Vl + Viz'+ 0{z^)] 
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with: 
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In typical appHcations, one is not interested in the physical distance D{z), but other definitions: 
• the luminosity distance: 
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the angular-diameter distance: 
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' sin(/;:_^ f(ff) k = +l; 
ro(i))= k = 0; 

. sinh(/;:_^ f(ff) fc = -l. 

If we make the expansion for short distances, namely if we insert the series expansion of a{t) in ro(-D), we have: 
dt 
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To convert from physical distance travelled to r coordinate traversed we have to consider that the Taylor series 
expansion of sin-sinh functions is: 
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so that Eq. lfTT]) with curvature k term becomes: 
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Using these one for luminosity distance we have: 
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While for the angular diameter distance it is: 
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K we want to use the same notation of [4011 , we define flo — 1 + -M-^ , which can be considered a purely cosmographic 

parameter, ot Qq — I — Qk — ^m.,o + ^r,o + ^x,o if we consider the dynamics of the Universe. With this parameter 
Eqs.(26)-(28) become: 

= l (53) 

+ (54) 



and 



Vl^^^-{l2-5qo + 3ql-Jo~no) (55) 
^i,a = ^ [52 - 20go + 21<7o' - 15<7^ - 7jo + lOgoio + - 2r!o(l + Sgo)] (56) 

2?i,^ = Y^ [359 - 184<7o + 186ql - ISbq^o + 105g^ + 90goJo - lOdq^jo - 15goso+ (57) 
- 57io + lOf + 9so -lo- 5r!o(17 - 690 + 9g§ - 2jo)] (58) 

25a,, = 1 (59) 
25Ay = -^(l + 9o) (60) 

'DXy = -l[-1o-3<ll+Ja + no] (61) 

2^12; = -^ [-2go + 3go' + 15g|5 - jo - Wqojo -so + 2^0] (62) 

'^Xy^-^ [1 - 6go + 9qo - 15q^o - ^^bq^ + 10 Wo + lOSg^jo + 15qoSo+ (63) 
- 3jo - 10 f + S0 + I0 + 50o] (64) 

Previous relations in this section have been derived for any value of the curvature parameter; but since in the 
following we will assume a fiat Universe, we will used the simplified versions for k = 0. Now, since we are going to 
use supernovae data, it will be useful to give as well the Taylor series of the expansion of the luminosity distance at it 
enters the modulus distance, which is the quantity about which those observational data inform. The final expression 
for the modulus distance based on the Hubble free luminosity distance, /i(z) = 5 log^o '^l(2), is: 

fi{z) ^ • (logz + M^z + M^z^ + M^z^ + M*z*) , (65) 



with 



M^ = -^[-l + qo] , (66) 

M^ = -j^[7- lOgo - 9q^o + 4io] , (67) 

M' = ^ [5 - 9go - 16ql - lOgg + 7jo + Sqojo + sq] , (68) 

M'^^-^ [-469 + 1004^0 + 2654g2 + 3300^3 + 1575^4 ^ 200^2 _ 1148jo+ 
2880 

- -2620goio - ISOOg^jo - 300goSo - 324so - 24^0] • (69) 



III. /(i?)-GRAVITY VS COSMOGRAPHY 
A. f{R) preliminaries 

As discussed in the Introduction, much interest has been recently devoted to the possibility that dark energy could 
be nothing else but a curvature effect according to which the present Universe is filled by pressureless dust matter 
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only and the acceleration is the result of modified Friedmann equations obtained by replacing the Ricci curvature 
scalar R with a generic function /(i?) in the gravity action. Under the assumption of a flat Universe, the Hubble 
parameter is therefore determined by^ : 



3 



f'iR) 



(70) 



where the prime denotes derivative with respect to R and Pcurv is the energy density of an effective curvature fl,uid^ : 



Pcurv = |i UiR) - Rf'm - 3H7?/"(i?)| . (71) 

Assuming there is no interaction between the matter and the curvature terms (we are in the so-called Jordan frame) , 
the matter continuity equation gives the usual scaling pi\f — pM{t = to)a^^ = iH^VLMa"^ , with VIm the present day 
matter density parameter. The continuity equation for pcurv then reads : 



611 [1 + Wcurv) Pcurv — : ; — — t « I'^j 



with 



= -1 



Rf'iR) + R Rf"'{R) - Hf"{R) 



[f{R)-Rf'iR)]/2-3HRf"iR) 



(73) 



the barotropic factor of the curvature fluid. It is worth noticing that the curvature fluid quantities pcurv and Wcurv 
only depends on f{R) and its derivatives up to the third order. As a consequence, considering only their present 
day values (which may be naively obtained by replacing R with Rq everywhere), two f{R) theories sharing the same 
values of f{Ro), f'{Ro), /"(i?o), /'"{Rq) will be degenerate from this point of view^. 

Combining Eq. (|72|) with Eq. (|70|) . one finally gets the following master equation for the Hubble parameter : 



H = -jj^{^H'o^Ma-^ + Rf"iR)- 



-R 



Rr"iR)-Hf"iR)]] 



(74) 



Expressing the scalar curvature R as function of the Hubble parameter as : 



i? -6 (i7 + 2H^^ (75) 

and inserting the result into Eq. l|74p . one ends with a fourth order nonlinear differential equation for the scale factor 
a{t) that cannot be easily solved also for the simplest cases (for instance, f{R) oc i?"). Moreover, although technically 
feasible, a numerical solution of Eq. (|74|) is plagued by the large uncertainties on the boundary conditions (i.e., the 
present day values of the scale factor and its derivatives up to the third order) that have to be set to find out the 
scale factor. 



^ We use here natural units such that S-rrG = 1. 

* Note that the name curvature fluid does not refer to the FRW curvature parameter k, but only takes into account that such a term is 
a geometrical one related to the scalar curvature R. 

^ One can argue that this is not strictly true since different f{R) theories will lead to different expansion rate H{t) and hence different 
present day values of R and its derivatives. However, it is likely that two f{R) functions that exactly match each other up to the third 
order derivative today will give rise to the same H{t) at least for t ~ to so that (-Ro, Ro, Ro) will be almost the same. 
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B. /(7?)-derivatives and cosmography 

Motivated by these difficulties, we approach now the problem from a different viewpoint. Rather than choosing a 
parameterized expression for f{R) and then numerically solving Eq. l(74|) for given values of the boundary conditions, we 
try to relate the present day values of its derivatives to the cosmographic parameters {qo,jo, sqJo) so that constraining 
them in a model independent way gives us a hint for what kind of f{R) theory could be able to fit the observed Hubble 
diagram®. 

As a preliminary step, it is worth considering again the constraint equation (|75l) . Differentiating with respect to t, 
we easily get the following relations : 

R = -6(h + iHH^ 

R = -6 (d^H/dt^ + AHH + 4ij2) . (76) 

d^R/dt^R = -6 {d^H/dt^ + AHd^H/dt^ + 12HH^ 
Evaluating these at the present time and using Eqs.® - Q, one finally gets : 

Rq ^ -&Hl{l - qo) , (77) 
i?o = -6ff3(jo - go - 2) , (78) 
/?o = -6iJo (so + <Zo + 8^0 + 6) , (79) 



d'^Ro/df' = ~QHl [lo - so + 2{qo + 4)jo - 6(3go + 8)go - 24] , (80) 

which will turn out to be useful in the following. 

Let us now come back to the expansion rate and master equations ((TOl) and l(74|) . Since they have to hold along 
the full evolutionary history of the Universe, they naively hold also at the present day. As a consequence, we may 
evaluate them in t = to thus easily obtaining : 

2 _ H^flM f{Ro) - RofiRp) - GHoRof"{Ro) 
°^/'(i?o)+ 6/'(i?o) ' 

SH^Qm ^o/"'(-Ro) + (^0 - HqRo^ /"{Rq) 

-^" = 2r(i?^+ wm • ^''^ 

Using Eqs.® - Q and ijTTll - 1(80|) . we can rearrange Eqs.JHT]) and (|82)) as two relations among the Hubble constant 
Ho and the cosmographic parameters (go, jo, sq), on one hand, and the present day values of f{R) and its derivatives 
up to third order. However, two further relations are needed in order to close the system and determine the four 
unknown quantities f{Ro), f'{Ro), /"(-Ro), f" {Ro)- A first one may be easily obtained by noting that, inserting back 
the physical units, the rate expansion equation reads : 



Note that a similar analysis, but in the context of the energy conditions in f(R), has yet been presented in |122l |. However, in that work, 
the author give an expression for f{R) and then compute the snap parameter to be compared to the observed one. On the contrary, 
our analysis does not depend on any assumed functional expression for f{R). 
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SttG 



[f-'rn + Pcurvf {R)] 



which clearly shows that, in f{R) gravity, the Newtonian gravitational constant G is replaced by an effective (time 
dependent) G^ff = G/f'{R). On the other hand, it is reasonable to assume that the present day value of Gg// is the 
same as the Newtonian one so that we get the simple constraint : 



Ge//(^ = 0) = G^/'(i?o) = l 



(83) 



In order to get the fourth relation we need to close the system, we first differentiate both sides of Eq. ((74l) with respect 
to t. We thus get : 



H = 



R^f"{R) + (i? - HR^ f"{R) + iH^flMa-^ 7?3/(»)(i?) + (3RR - HR^^ f"{R) 



R!"{R) [.f'{R)Y 



2/'(i?) 



SR/dt^ - HR + HRj f"{R) - m^flMHa 



(84) 



with /^™^ (R) = d'^f/dR'^. Let us now suppose that /(i?) may be well approximated by its third order Taylor expansion 
in R — Rq, i.e. we set : 



f{R) - f{Ro) + /'(i?o)(i? - Ro) + \f"{R^){R - Rof + \f"'{Ro){R - Rof 



(85) 



In such an approximation, it is /("^(i?) — d^^f/R^^ = for n > 4 so that naively /'-"•'(-Ro) = 0. Evaluating then 
Eq. I|84p at the present day, we get : 



Ho 



Rlf"{Ra) + {Ro - HoRo) f"{Ro) + iHl^M (iiloRo - Hktj f"{Ro) 



Ra.f"{Ra) [r{Ro)f 



2/'(i?o) 



(^d^Ro/dt^ - HoRo + HoRo) f"{Ro) - ^H^^m 



V'iRo) 



(86) 



We can now schematically proceed as follows. Evaluate Eqs.® - Q at 2 = and plug these relations into the left 
hand sides of Eqs. ljST]) . l(82|). (|86l) . Insert Eqs. l(77|) - (|80|) into the right hand sides of these same equations so that 
only the cosmographic parameters {qo, ja, so,la) and the f{R) related quantities enter both sides of these relations. 
Finally, solve them under the constraint l(83l) with respect to the present day values of f{R) and its derivatives up to 
the third order. After some algebra, one ends up with the desired result : 



fjRo) 
ml 



'Po{qo,jo,soJo)^M + Qoiqo,jo,so,lo) 
T^iqoJo^so^k) 



(87) 



f{Ro) = 1 , 



f"{Ro) _ 'P2{qo,jo,So)i^M + £2(90, jo, gp) 

[GH^y^ 7^((7o,jo,so,^o) 



(89) 



f"'{Ro) _ V'i{qo,io,so,lo)^M + Qi{qo,io,so.,lo) 



mi] 



(jo - <?o - 2)7^(go, jo, So, ^o) 



(90) 
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where we have defined : 

= Uo - 90 ~ 2)^0 - (3so + 7ja + Gql + 41<7o + 22)sq - [{Sqo + 16)jo + 20qi + 64go + 12] jo + 



- (3g;5 + 25%' + 96^0 + 72go + 20) , (91) 

Qo = {ql - mo + 2go)^o + [3qoSo + {Mo + 6)jo + Qqo + 44^^ + 22go - 12] sq 

+ [2j2 + (3g2 + lOgo - 6)io + I7ql + 52ql + 54qo + 36] jo + 3ql + 28ql + U8q^ + 

+ 72ql - 76<7o - 64 , (92) 

V2 = 9so + 6io + %l + 66go + 42 , (93) 

22 = -{6(go + l)so + [2jo-2(l-go)]jo + 6go' + 50g2 + 749o + 32} , (94) 

P3 = 3/o + 3so - 9(go + 4)io - (45g2 + 78qo + 12) , (95) 

23 - - {2(1 + go)^o + 2(jo + 9o)so - (2jo + 4go' + 12go + 6) jo - (30g3 + Mql + 78go + 24)} (96) 

= Uo -qo- 2)/o - (3so - 2 jo + 6ql + 5Qqo + 40)so + [{3qo + 10) jo + 11^0 + 4^0+ 

- 18] jo - i'Sqt + 34g3 + 246go + 104) . (97) 



Eqs. (|87l) - 197|) make it possible to estimate the present day values of f{R) and its first three derivatives as function 
of the Hubble constant Hq and the cosmographic parameters (qo, jo, so,lo) provided a value for the matter density 
parameter VIm is given. This is a somewhat problematic point. Indeed, while the cosmographic parameters may be 
estimated in a model independent way, the fiducial value for flM is usually the outcome of fitting a given dataset in the 
framework of an assumed dark energy scenario. However, it is worth noting that different models all converge towards 
the concordance value flM — 0.25 which is also in agreement with astrophysical (model independent) estimates from 
the gas mass fraction in galaxy clusters. On the other hand, it has been proposed that /(i?) theories may avoid the 
need for dark matter in galaxies and galaxy clusters [iB, [H, HI, [13, 111, llOSl . Il46l |. In such a case, the total matter 
content of the Universe is essentially equal to the baryonic one. According to the primordial elements abundance and 
the standard BBN scenario, we therefore get flM — ^b/h'^ with Ui, = Vliih? ~ 0.0214 [93] and h the Hubble constant 
in units of lOOkm/s/Mpc. Setting h = 0.72 in agreement with the results of the HST Key project [t^I, we thus get 
r^M = 0.041 for a baryons only Universe. We will therefore consider in the following both cases when numerical 
estimates are needed. 

It is worth noticing that Ho only plays the role of a scaling parameter giving the correct physical dimensions to f{R) 
and its derivatives. As such, it is not surprising that we need four cosmographic parameters, namely (go, Jo, sq, ^o), 
to fix the four f{R) related quantities /(i?o), /'(^o), f"{Ro), f"'{Ro)- It is also worth stressing that Eqs. ([57)) - IpOl) 
are linear in the f{R) quantities so that (go, Jo, so, ^o) uniquely determine the former ones. On the contrary, inverting 
them to get the cosmographic parameters as function of the f{R) ones, we do not get linear relations. Indeed, the 
field equations in f{R) theories are nonlinear fourth order differential equations in the scale factor a(t) so that fixing 
the derivatives of f{R) up to third order makes it possible to find out a class of solutions, not a single one. Each 
one of these solutions will be characterized by a different set of cosmographic parameters thus explaining why the 
inversion of Eqs. i|87p - i(97|) does not give a unique result for (go, jo, so, ^o)- 

As a final comment, we reconsider the underlying assumptions leading to the above derived relations. While Eqs. (|8T|) 
and (|82l) are exact relations deriving from a rigorous application of the field equations, Eq. (|86|) heavily relies on having 
approximated f{R) with its third order Taylor expansion l(85|) . If this assumption fails, the system should not be 
closed since a fifth unknown parameter enters the game, namely /^™''(i?o)- Actually, replacing f{R) with its Taylor 
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expansion is not possible for all class of f{R) theories. As such, the above results only hold in those cases where such 
an expansion is possible. Moreover, by truncating the expansion to the third order, we are impHcitly assuming that 
higher order terms are negligible over the redshift range probed by the data. That is to say, we are assuming that : 



/(")(i?o)(i? - i?o)" « V ^ V^(i? - Ror for n > 4 (98) 

■^-^ m! 

m— 

over the redshift range probed by the data. Checking the validity of this assumption is not possible without explicitly 
solving the field equations, but we can guess an order of magnitude estimate considering that, for all viable models, 
the background dynamics should not differ too much from the ACDM one at least up to z ~ 2. Using then the 
expression of H{z) for the ACDM model, it is easily to see that R/Rq is a quickly increasing function of the redshift 
so that, in order Eq. ljMl) holds, we have to assume that /^"^(i?o) << f"'{Ro) for n > A. This condition is easier to 
check for many analytical f{R) models. 

Once such a relation is verified, we have still to worry about Eq. l|83p relying on the assumption that the cosmological 
gravitational constant is exactly the same as the local one. Although reasonable, this requirement is not absolutely 
demonstrated. Actually, the numerical value usually adopted for the Newton constant G^v is obtained from laboratory 
experiments in settings that can hardly be considered homogenous and isotropic. As such, the spacetime metric in 
such conditions has nothing to do with the cosmological one so that matching the two values of G is strictly speaking 
an extrapolation. Although commonly accepted and quite reasonable, the condition Giocai — Gcosmo could (at least, 
in principle) be violated so that Eq. ljM)) could be reconsidered. Indeed, as we will see, the condition /'(i?o) = 1 may 
not be verified for some popular f{R) models recently proposed in literature. However, it is reasonable to assume 
that Geff{z = 0) = G(l + e) with e << 1. When this be the case, we should repeat the derivation of Eqs. (|87| - 1(90|) 
now using the condition f'{Ro) = (1 + e)^^. Taylor expanding the results in e to the first order and comparing with 
the above derived equations, we can estimate the error induced by our assumption £ = 0. The resulting expressions 
are too lengthy to be reported and depend in a complicated way on the values of the matter density parameter Qm , 
the cosmographic parameters (go, jo, so, ^o) and e. However, we have numerically checked that the error induced on 
/(i?o), /"(i?o), f"'{Ro) are much lower than 10% for value of e as high as an unrealistic e ^ 0.1. We are confident 
that our results are reliable also for these cases. 



IV. /(i?)-GRAVITY AND THE CPL MODEL 

A determination of f{R) and its derivatives in terms of the cosmographic parameters need for an estimate of these 
latter from the data in a model independent way. Unfortunately, even in the nowadays era of precision cosmology, such 
a program is still too ambitious to give useful constraints on the f{R) derivatives, as we will see later. On the other 
hand, the cosmographic parameters may also be expressed in terms of the dark energy density and EoS parameters so 
that we can work out what are the present day values of f{R) and its derivatives giving the same (go, jo, so, ^o) of the 
given dark energy model. To this aim, it is convenient to adopt a parameterized expression for the dark energy EoS 
in order to reduce the dependence of the results on any underlying theoretical scenario. Following the prescription 
of the Dark Energy Task Force j2], we will use the Chevallier - Polarski - Linder (CPL) parameterization for the EoS 
setting (H, iil : 

W = Wo + Wa(l — a) = Wq + WaZ{l + z)~^ (99) 

so that, in a fiat Universe filled by dust matter and dark energy, the dimensionless Hubble parameter E{z) = H/Ha 
reads : 



E^{z) = nM{l + zf + nx{l + z)3(i+-o+-a)e-^ (100) 

with fix — 1 — i^Af because of the fiatness assumption. In order to determine the cosmographic parameters for such 
a model, we avoid integrating H{z) to get a{t) by noting that d/dt — —(1 + z)H{z)d/dz. We can use such a relation 
to evaluate {H, H, d^H/dt^, d'^H/dt'^) and then solve Eqs.lHl) - ([9]), evaluated in z = 0, with respect to the parameters 
of interest. Some algebra finally gives : 



qo^i^ + i^{l-^M)wo , (101) 
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jo = 1 + 2 (1 ^ ^Af) [3wo(l + Wo) + Wa] , 



(102) 



So 



33 



27, 



(1 - nM)Wa - -(1 - ^m) [9 + (7 - nM)Wa] Wo-'il- nM){l6 ~ 31}Af)u)2 - —(1 - Qm^S - ^M^ij^] 



35 

y 



[213 + (7 - ilM)Wa\ Wa H 



67-2117 



M 



-(23- I^mK 



81 



(i-r!M)(3-2r!MK 



[489 + 9(82 - 2in M)wa\ wq + 
27 

f — (l-l]M)(47- 24r!M)w'o 



(104) 



Inserting Eqs. ljlOip - l|104p into Eqs. l(87|) - (|97| . we get lengthy expressions (which we do not report here) giving the 
present day values of f{R) and its first three derivatives as function of {Qm ,WQ,Wa)- It is worth noting that the f{R) 
model thus obtained is not dynamically equivalent to the starting CPL one. Indeed, the two models have the same 
cosmographic parameters only today. As such, for instance, the scale factor is the same between the two theories only 
over the time period during which the fifth order Taylor expansion is a good approximation of the actual a{t). It is 
also worth stressing that such a procedure does not select a unique f{R) model, but rather a class of fourth order 
theories all sharing the same third order Taylor expansion of f{R). 



A. The ACDM case 

With these caveats in mind, it is worth considering first the ACDM model which is obtained by setting {wq, Wa) 
(—1,0) in the above expressions thus giving: 



1 3 

90 = 2^2*^^^ ^^''^ 



jo = 1 



So — I — -Hm 

27 

lo = 1 + 3r2j\/ + — rijyj 



(105) 



When inserted into the expressions for the f{R) quantities, these relations give the remarkable result : 



f{Ro)=Ro + 2A , /"(i?o) = /"'(i?o) = 



(106) 



so that we obviously conclude that the only f{R) theory having exactly the same cosmographic parameters as the 
ACDM model is just f{R) oc R, i.e. GR. It is worth noticing that such a result comes out as a consequence of the 
values of (go, Jo) in the ACDM model. Indeed, should we have left (sq, Iq) undetermined and only fixed {qo,jo) to the 
values in l|105p . we should have got the same result in (|106p . Since the ACDM model fits well a large set of different 
data, we do expect that the actual values of {qo, jo, so,lo) do not differ too much from the ACDM ones. Therefore, 
we plug into Eqs. fSTj) - (|97| the following expressions : 



qo = qax{l + eq) , jo=j^x(l + £j) 



so = sgx{l + es) , lo^lo^il + ei) 
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with {qq , Jq , ,Iq) given by Eqs. l|105p and (e^, e^, Eg, £;) quantifyin the deviations from the ACDM values allowed 
by the data. A numerical estimate of these quantities may be obtained, e.g., from a Markov chain analysis, but 
this is outside our aims. Since we are here interested in a theoretical examination, we prefer to consider an ideaHzed 
situation where the four quantities above all share the same value e << 1. In such a case, we can easily investigate how 
much the corresponding f{R) deviates from the GR one considering the two ratios /"(i?o)//(-Ro) and f'"{Ro)/f{Ro)- 
Inserting the above expressions for the cosmographic parameters into the exact (not reported) formulae for /(i?o), 
f"{Ro) and /"'(i?o), taking their ratios and then expanding to first order in e, we finally get : 



64-6»m(9»m + 8) £_ 
[3(90m + 74)nM - 556] + 16 ^ 27 ' 



_ 6 [{smM - iio)»M + 40] nM + 16 £ 

'^^^ [3{9nM + 7A)nM - 556] ftj^j + 16 243172^ ' ^ ' 

having defined 7720 = f"{Ro)/f{Ro)xHQ and 7730 — f"'{Ro)/f{Ro)xHQ which, being dimensionless quantities, are 
more suited to estimate the order of magnitudes of the different terms. Inserting our fiducial values for ^Im, we get : 

7720 ^ 0.15 X £ for Qm = 0.041 

7^20 ^ -0.12 X £ for Qm = 0.250 

7730 ~ 4 X £ for Qm = 0.041 
7730 ~ -0.18 X £ for Qm = 0.250 

For values of £ up to 0.1, the above relations show that the second and third derivatives are at most two orders of 
magnitude smaller than the zeroth order term /(i?o)- Actually, the values of 7^30 for a baryon only model (first row) 
seems to argue in favor of a larger importance of the third order term. However, we have numerically checked that 
the above relations approximates very well the exact expressions up to £ ~ 0.1 with an accuracy depending on the 
value of flM, being smaller for smaller matter density parameters. Using the exact expressions for 7720 and 7730, our 
conclusion on the neghgible effect of the second and third order derivatives are significantly strengthened. 

Such a result holds under the hypotheses that the narrower are the constraints on the validity of the ACDM model, 
the smaller are the deviations of the cosmographic parameters from the ACDM ones. It is possible to show that this 
indeed the case for the CPL parametrization we are considering. On the other hand, we have also assumed that the 
deviations (£q, £j, £s, ei) take the same values. Although such hypothesis is somewhat ad hoc, we argue that the main 
results are not affected by giving it away. Indeed, although different from each other, we can still assume that all of 
them are very small so that Taylor expanding to the first order should lead to additional terms into Eqs. ljlOTp - IjlOSp 
which are likely of the same order of magnitude. We may therefore conclude that, if the observations confirm that 
the values of the cosmographic parameters agree within ~ 10% with those predicted for the ACDM model, we must 
conclude that the deviations of f{R) from the GR case, f{R) cx R, should be vanishingly small. 

It is worth stressing, however, that such a conclusion only holds for those /(i?) models satisfying the constraint 
((98| . It is indeed possible to work out a model having f{Ro) oc i?o, f"{Ro) = f'"{Ro) = 0, but /'^"^(i?o) ^ for some 
n. For such a (somewhat ad hoc) model, Eq.ljOS]) is clearly not satisfied so that the cosmographic parameters have to 
be evaluated from the solution of the field equations. For such a model, the conclusion above does not hold so that 
one cannot exclude that the resulting {qo, jo, sqJq) are within 10% of the ACDM ones. 



B. The constant EoS model 



Let us now take into account the condition w = —1, but still retains Wa = thus obtaining the so called quiessence 
models. In such a case, some problems arise because both the terms (jo — go — 2) and TZ may vanish for some 
combinations of the two model parameters {flM,wo)- For instance, we find that jo — 90 — 2 = for wq = {wi,W2) 
with : 

1 

^1 = = , 

1 - + V(i-^^M)(4-r!M) 



15 



0.01 
0.005 


= -0.005 
-0.01 
-0.015 
-0.02 



-1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 

Wo 

Figure 1: The dimensionless ratio between the present day values of f"[R) and f[R) as function of the constant EoS luo of the 
corresponding quiessence model. Short dashed and solid lines refer to models with Q,m = 0.041 and 0.250 respectively. 



4- »Af 

v/(i-f^M)(4-r!M)_ 

On the other hand, the equation TZ{flM, ^o) — may have different real roots for w depending on the adopted value 
of SIm- Denoting collectively with w„„h the values of wo that, for a given Qm, make (jo — qo — 2)7^(riA/, wq) taking 
the null value, we individuate a set of quiessence models whose cosmographic parameters give rise to divergent values 
of f{Ro, f"{Ro) and /"'(i?o)- For such models, f{R) is clearly not defined so that we have to exclude these cases 
from further consideration. We only note that it is still possible to work out a f{R) theory reproducing the same 
background dynamics of such models, but a different route has to be used. 

Since both go and jo now deviate from the ACDM values, it is not surprising that both /"(i?o) and f"'{Ro) take 
finite non null values. However, it is more interesting to study the two quantities 7720 and 7730 defined above to 
investigate the deviations of f{R) from the GR case. These are plotted in Figs. [T] and [2] for the two fiducial SI a/ values. 
Note that the range of wo in these plots have been chosen in order to avoid divergences, but the lessons we will draw 
also hold for the other wq values. 

As a general comment, it is clear that, even in this case, /"(i?o) and /"'(i?o) are from two to three orders of 
magnitude smaller that the zeroth order term /(i?o)- Such a result could be yet guessed from the previous discussion 
for the ACDM case. Actually, relaxing the hypothesis wo = — 1 is the same as allowing the cosmographic parameters 
to deviate from the ACDM values. Although a direct mapping between the two cases cannot be estabHshed, it is 
nonetheless evident that such a relation can be argued thus making the outcome of the above plots not fully surprising. 
It is nevertheless worth noting that, while in the ACDM case, 7720 and 7730 always have opposite signs, this is not the 
case for quiessence models with w > —1. Indeed, depending on the value of f^A/, we can have f{R) theories with both 
7720 and 7730 positive. Moreover, the lower is Qm, the higher are the ratios 7720 and 7730 for a given value of wq. This 
can be explained qualitatively noticing that, for a lower fin, the density parameter of the curvature fiuid (playing 
the role of an effective dark energy) must be larger thus claiming for higher values of the second and third derivatives 
(see also [28| for a different approach to the problem) . 

C. The general case 

Finally, we consider evolving dark energy models with Wa 9^ 0. Needless to say, varying three parameters allows 
to get a wide range of models that cannot be discussed in detail. Therefore, we only concentrate on evolving dark 
energy models with wq = — 1 in agreement with some most recent analysis. The results on 7720 and 7730 are plotted in 
Figs.[3]and[4]where these quantities as functions of Wa- Note that we are considering models with positive Wa so that 
w{z) tends to wq + Wa > wq for z ^ 00 so that the EoS dark energy can eventually approach the dust value w = 0. 
Actually, this is also the range favored by the data. We have, however, excluded values where 7720 or 7730 diverge. 
Considering how they are defined, it is clear that these two quantities diverge when /(i?o) = so that the values of 
{wojWa) making (7720, ^o) to diverge may be found solving: 




1 

^2 = -3 



Vo{wo,Wa)flM + Qo{'Wo,Wa) = 
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Figure 2: The dimensionless ratio between the present day values of f"'{R) and f{R) as function of the constant EoS wo of 
the corresponding quiessence model. Short dashed and solid lines refer to models with Q,m = 0.041 and 0.250 respectively. 
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Figure 3: The dimensionless ratio between the present day values of f"{R) and f{R) as function of the Wa parameter for models 
with Wo — —1. Short dashed and solid lines refer to models with ^Im = 0.041 and 0.250 respectively. 

where VQ{wQ,Wa) and QQ{wo,Wa) are obtained by inserting Eqs. (|10ip - l|104p into the defintions (|9T|) - 192|) . For such 
CPL models, there is no any f{R) model having the same cosmographic parameters and, at the same time, satisfying 
all the criteria needed for the validity of our procedure. Actually, if f{Ro) — 0, the condition (|98|) is likely to be 
violated so that higher than third order must be included in the Taylor expansion of f{R) thus invalidating the 
derivation of Eqs.(l87]) - ([90)) . 

Under these caveats. Figs. [3] and |4] demonstrate that alloviring the dark energy EoS to evolve does not change 
significantly our conclusions. Indeed, the second and third derivatives, although being not null, are nevertheless 
negligible with respect to the zeroth order term thus arguing in favour of a GR-like f{R) with only very small 
corrections. Such a result is, however, not fully unexpected. From Eqs. (|10ip and (|102p . we see that, having setted 
Wo = —1, the go parameter is the same as for the ACDM model, while jo reads Jq + (3/2)(l — il,M)wa. As we have 
stressed above, the Hilbert - Einstein Lagrangian f{R) = i? + 2A is recovered when {qo,jo) — (QajJo) whatever 
the values of {sqJq) are. Introducing a Wa ^ makes {sqJo) to differ from the ACDM values, but the first two 
cosmographic parameters are only mildly affected. Such deviations are then partially washed out by the complicated 
way they enter in the determination of the present day values of f{R) and its first three derivatives. 

V. CONSTRAINING f{R) PARAMETERS 

In the previous section, we have worked an alternative method to estimate /(i?o), f'iRo), f"'{Ra) resorting to a 
model independent parameterization of the dark energy EoS. However, in the ideal case, the cosmographic parameters 
are directly estimated from the data so that Eqs. l(87|) - ((97l) can be used to infer the values of the f{R) related quantities. 
These latter can then be used to put constraints on the parameters entering an assumed fourth order theory assigned 
by a f{R) function characterized by a set of parameters p = (pi, . . . ,p„) provided that the hypotheses underlying 
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Figure 4: The dimensionless ratio between the present day values of f'"{R) and fiR) as function of the Wa parameter for 
models with wo = —1. Short dashed and solid lines refer to models with = 0.041 and 0.250 respectively. 



the derivation of Eqs. ((87|l - (|97)) are indeed satisfied. We show below two interesting cases which clearly highlight the 
potentiality and the limitations of such an analysis. 



A. Double power law Lagrangian 



As a first interesting example, we set : 



f{R) =R{1 + ai?" + (109) 

with n and m two positive real numbers (see, for example, [ll6l | for some physical motivations). The following 
expressions are immediately obtained : 

' /(i?o) = Ro {1 + aB^ + fSRo"") 

f{Ro) = 1 + a(n + l)i?ff - /3(m - l)i?,o " 

' f'{Ro) = an{n + l)R^-^ +Pm{m-1)R^'^^+"'^ ' 

f"{Ro) = an(n+l)(n-l)i?^2 

- /3m(m + l)(m-l)i?o('+"^ 

Denoting by 0; (with i = 0, . . . , 3) the values of (Ro) determined through Eqs. (|87l) - (|97| . we can solve : 

r f{Ro) = <f>o 

^ /'(i?o) = 01 

/"(i?0) = 02 
[ f'iRo) = 03 

which is a system of four equations in the four unknowns {a,(3,n,m) that can be analytically solved proceeding as 
follows. First, we solve the first and second equation with respect to {a, (3) obtaining: 



1 — m 
n + m 

l + n 



00 

Ro 



Rn 



(110) 



n + m\ i?o ' 



while, solving the third and fourth equations, we get : 
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a 



f3 



cl)2Ro[l+m + i4>3/^2)Ro] 
n(n + l)(n + m) 



hRl+''[l-n+{cj)3/^2)Ro] 



(111) 



m(l — m){n + m) 

Equating the two solutions, we get a systems of two equations in the two unknowns (n, m), namely : 



n(n+l)(l-m)(l-(/)o/flo) 
cj)2Ro [l + m+(03/02)-Ro] 

m(n + l)(m-l)(l-0o/i?o) 



= 1 



(112) 



iRo [1 - n 



I^3/<P2)Rq 



= 1 



Solving with respect to m, we get two solutions, the first one being m 
(a, P) goes to infinity. The only acceptable solution is : 



-n which has to be discarded since makes 



m = -[l-n+(03/02)i?o] (113) 
which, inserted back into the above system, leads to a second order polynomial equation for n with solutions : 



where we have defined : 



1 



i + — /to±- 



2Ro{1 + (I>q/Ro) 



(114) 



A/'((/)o>2,03) = 



{Rlcbl - 2i?3</.o + Rt) cbj 
6 (i?o0o - 2i?§0o + Rl) <l)2h 
9 (</.2 - 27?o(/.o + Rl) 4>l 
4 (Rlcj^o - Rl) 4>l . 



(115) 



Depending on the values of (go, Jo, Sq, /o), Eq. l|114p may lead to one, two or any acceptable solution, i.e. real positive 
values of n. This solution has then to be inserted back into Eq. (|113p to determine m and then into Eqs. l|110p or l|llip 
to estimate (a, (3). If the final values of (a, (}, n, m) are physically viable, we can conclude that the model in Eq. l|109p 
is in agreement with the data giving the same cosmographic parameters inferred from the data themselves. Exploring 
analytically what is the region of the (go, io, so, ^o) parameter space which leads to acceptable (a, /3, n, m) solutions is 
a daunting task far outside the aim of the present work. 



B. The Hu and Sawicki model 



One of the most pressing problems of f{R) theories is the need to escape the severe constraints imposed by the 
Solar System tests. A successful model has been recently proposed by Hu and Sawicki [s^ (HS) setting^ : 



f{R) = R-Ra 



ajR/RcT 
1 + P{R/RcT 



(116) 



As for the double power law model discussed above, there are four parameters which we can be expressed in terms of 
the cosmographic parameters (go, io, so, ^o)- 



Note that such a model does not pass the matter instabiUty test so that some viable generalizations [i?! . Ill4l . lll5l | have been proposed. 
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As a first step, it is trivial to get : 



with Roc — Rq/Rc and : 



fiRo) = 
f'{Ro) = 

/"(i?o) = 

f"{Ro) = 
A 



Rq — R. 



1 + m'^c 

anRcRpc 
' R^{l + l3Roc? 

anRcRSc [(1 - n) + + n)R^^ 

RU^ + m^c)' 

anRcR}i^{An^ + Bn + C) 

rW+MJ 



(117) 



/3^i?^^ + 4/3i?^,-l 



B = 3(1 



P'Rlc) 



(118) 



c 



-2(1 - PR^, 



Equating Eqs. (|117p to the four quantities {4>Q,(l)i,4>2,4>'i) defined as above, we could, in principle, solve this system 
of four equations in four unknowns to get {a, P, Rc,n) in terms of {(f)o , , (j)2 , (j)3) and then, using Eqs. (|87| - Ip7|l as 
functions of the cosmographic parameters. However, setting = 1 as required by Eq. l|88p gives the only trivial 
solution anRc = so that the HS model reduces to the Einstein - Hilbert Lagrangian f{R) — R. In order to escape 
this problem, we can relax the condition f'{Ro) = 1 to f'{Ro) = (1 + s)^^ ■ As we have discussed in Sect. IV, this is 
the same as assuming that the present day effective gravitational constant Ge//,o — Gn / f [Ro) only slightly differs 
from the usual Newtonian one which seems to be a quite reasonable assumption. Under this hypothesis, we can 
analytically solve for (a, (3, Rc, n) in terms of (0o, £, </'2, ^a)- The actual values of (00, 02, ^a) will be no more given by 
Eqs. (|87l) - 1(90|) . but we have checked that they deviate from those expressions* much less than 10% for e up to 10% 
well below any realistic expectation. 
With this caveat in mind, we first solve 



/(i?o) = 0o , /"(i?o) = (1 + e)" 



to get : 



/3 



n(l + £) /i?o 
e \Rc 
n{l+e) (Rq 
Rc 



1 _ 



i?0 



1_ 



i?o n(l + e) 



Inserting these expressions in Eqs. (|117p . it is easy to check that Rc cancels out so that we can no more determine its 
value. Such a result is, however, not unexpected. Indeed, Eq. l|116p can trivially be rewritten as: 



f{R) = R 



aR" 



1 + /3i?" 



with a = aRl " and /3 = f3R^ " which are indeed the quantities that are determined by the above expressions 
for (a,/3). Reversing the discussion, the present day values of /''•'(i?) depend on {a,P,Rc) only through the two 



Note that the correct expressions for (p/i-io , 02 , 03 ) may still formally be written as Eqs.JSTj - JSoJ, but the polynomials entering them 
are now different and also depend on powers of e. 
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parameters (a,/3). As such, the use of cosmographic parameters is unable to break this degeneracy. However, since 
Rc only plays the role of a scaling parameter, we can arbitrarily set its value without loss of generality. 

On the other hand, this degeneracy allows us to get a consistency relation to immediately check whether the HS 
model is viable or not. Indeed, solving the equation /"(i?o) = </'2, we get : 

(0o/j?o) + [(1 + £)/£](! - 02i^o) - (1 - £)/(! + e) 
1 - 0o/i?o 

which can then be inserted into the equations f"'{Ro) = 03 to obtain a complicated relation among {(j>o, (t>2,<t>-i) which 
we do not report for sake of shortness. Solving such a relation with respect to 03 /0o and Taylor expanding to first 
order in e, the constraint we get reads : 



1 



Rq 



00 



1 - 



2e 



1 - 0o/i?o 



If the cosmographic parameters (go, io: soj ^o) are known with sufficient accuracy, one could compute the values of 
(i?07 00: 02-03) for a given e (eventually using the expressions obtained for e = 0) and then check if they satisfied this 
relation. If this is not the case, one can immediately give off the HS model also without the need of solving the field 
equations and fitting the data. Actually, given the still large errors on the cosmographic parameters, such a test only 
remains in the realm of (quite distant) future applications. However, the HS model works for other tests as shown in 
[88] and so a consistent cosmography analysis has to be combined with them. 



VI. CONSTRAINTS ON /(i?)-DERIVATIVES FROM THE DATA 

Eqs. ljSTj) - (|97l) relate the present day values of f{R) and its first three derivatives to the cosmographic parameters 
(qo, jo, so,lo) and the matter density flM- In principle, therefore, a measurement of these latter quantities makes it 
possible to put constraints on f^^\Ro), with i = {0, ... ,3}, and hence on the parameters of a given fourth order 
theory through the method shown in the previous section. Actually, the cosmographic parameters are affected by 
errors which obviously propagate onto the f{R) quantities. Actually, the covariance matrix for the cosmographic 
parameters is not diagonal so that one has also take care of this to estimate the final errors on /'•'^(i?o)- A similar 
discussion also holds for the errors on the dimensionless ratios 7720 and 7730 introduced above. As a general rule, 
indicating with g{Q]\i,p) a generic f{R) related quantity depending on Qm and the set of cosmographic parameters 
p, its uncertainty reads : 



-1 = 



dg 



on 



M 



' M 



E 

i=l 



dg 
dpi 



dpi dpj 



(119) 



where are the elements of the covariance matrix (being Cu — CTp.), we have set {pi,P2,P3,P4) — (go, Jo, so, ^o)- and 
assumed that the error ctm on is uncorrelated with those on p. Note that this latter assumption strictly holds if 
the matter density parameter is estimated from an astrophysical method (such as estimating the total matter in the 
Universe from the estimated halo mass function). Alternatively, we will assume that 17m is constrained by the CMBR 
related experiments. Since these latter mainly probes the very high redshift Universe {z ~ ziss — 1089), while the 
cosmographic parameters are concerned with the present day cosmo, one can argue that the determination of is 
not affected by the details of the model adopted for describing the late Universe. Indeed, we can reasonably assume 
that, whatever is the dark energy candidate or /(i?) theory, the CMBR era is well approximated by the standard GR 
with a model comprising only dust matter. As such, we will make the simplifying (but well motivated) assumption 
that tJM niay be reduced to very small values and is uncorrelated with the cosmographic parameters. 

Under this assumption, the problem of estimating the errors on g{flM,p) reduces to estimating the covariance 
matrix for the cosmographic parameters given the details of the data set u sed as observational constraints. We 
address this issue by computing the Fisher information matrix (see, e.g., [l58{ and references therein) defined as: 



with L = —2lnC{9i, . . . ,9n), C{9i, . . . ,9n) the likeHhood of the experiment, (6'i,...,0„) the set of parameters to 
be constrained, and (. . .) denotes the expectation value. Actually, the expectation value is computed by evaluating 
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the Fisher matrix elements for fiducial values of the model parameters (6'i, . . . , 0„), while the covariance matrix C is 
finally obtained as the inverse of F. 

A key ingredient in the computation of F is the definition of the likelihood which depends, of course, of what 
experimental constraint one is using. To this aim, it is worth remembering that our analysis is based on fifth order 
Taylor expansion of the scale factor a{t) so that we can only rely on observational tests probing quantities that are well 
described by this truncated series. Moreover, since we do not assume any particular model, we can only characterize 
the background evolution of the Universe, but not its dynamics which, being related to the evolution of perturbations, 
unavoidably need the specification of a physical model. As a result, the SNela Hubble diagram is the ideal test^ to 
constrain the cosmographic parameters. We therefore defined the likelihood as : 



C{Ho,p)<xeiip-xHHo..p)/2 



where the distance modulus to redshift z reads : 



^J■obsiZi) - ^Jith{Zn, Ho, p) 



(121) 



fith{z,Ho,p) = 25 + 51og(c/i/o) + 51ogdL(z,p) 
and dL{z) is the Hubble free luminosity distance : 



(122) 



Using the fifth order Taylor expansion of the scale factor, we get for ^^(z, p) an analytical expression (reported in 
Appendix A) so that the computation of Fij does not need any numerical integration (which makes the estimate 
faster). As a last ingredient, we need to specify the details of the SNela survey giving the redshift distribution of the 
sample and the error on each measurement. Following f9^, we adopt^° : 




with Zmax the maximum redshift of the survey, asys an irreducible scatter in the SNela distance modulus and am to 
be assigned depending on the photometric accuracy. 

In order to run the Fisher matrix calculation, we have to set a fiducial model which we set according to the ACDM 
predictions for the cosmographic parameters. For = 0.3 and h = 0.72 (with h the Hubble constant in units of 
lOOkm/s/Mpc), we get : 

(go, Jo, so, ^o) = (-0.55,1.0,-0.35,3.11) . 

As a first consistency check, we compute the Fisher matrix for a survey mimicking the recent database in [g^I thus 
setting {NsNeia,<ym) — (192,0.33). After marginalizing over h (which, as well known, is fully degenerate with the 
SNela absolute magnitude 7W), we get for the uncertainties : 

(cTi, (72, (73, (T4) = (0.38, 5.4, 28.1, 74.0) 

where we are still using the indexing introduced above for the cosmographic parameters. These values compare 
reasonably well with those obtained from a cosmographic fitting of the Gold SNela dataset^^ (s^ . [83| : 



^ See the conclusions for further discussion on this issue. 
Note that, in [93 |. the authors assume the data are separated in redshift bins so that the error becomes tr-^ = ^lys/Mhin + 
J^bin(z / ^max)'^ o'^ with Nbin the number of SNela in a bin. However, we prefer to not bin the data so that Nun = 1- 
Actually, such estimates have been obtained computing the mean and the standard deviation from the marginalized likelihoods of the 
cosmographic parameters. As such, the central values do not represent exactly the best fit model, while the standard deviations do not 
give a rigorous description of the error because the marginalized likelihoods are manifestly non - Gaussian. Nevertheless, we are mainly 
interested in an order of magnitude estimate so that we do not care about such statistical details. 
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go = -0.90±0.65 , jo = 2.7±6.7 , 

So = 36.5±52.9 , ^o = 142.7±320 . 

Because of the Gaussian assumptions it relies on, the Fisher matrix forecasts are known to be lower Hmits to the 
accuracy a given experiment can attain on the determination of a set of parameters. This is indeed the case with the 
comparison suggesting that our predictions are quite optimistic. It is worth stressing, however, that the analysis in 
[82, 83] used the Gold SNela dataset which is poorer in high redshift SNela than the jfi^l one we are mimicking so 
that larger errors on the higher order parameters (sqj ^o) are expected. 

Rather than computing the errors on f{Ro) and its first three derivatives, it is more interesting to look at the 
precision attainable on the dimensionless ratios (7720, introduced above since they quantify how much deviations 
from the linear order are present. For the fiducial model we are considering, both 7720 and 7730 vanish, while, using the 
covariance matrix for a present day survey and setting om I^m — 10%, their uncertainties read : 

(^T20,fT3o) = (0.04,0.04) . 

As an application, we can look at Figs. [1] and [2] showing how (77207 ?73o) depend on the present day EoS wq for f{R) 
models sharing the same cosmographic parameters of a dark energy model with constant EoS. As it is clear, also 
considering only the lu range, the full region plotted is allowed by such large constraints on (7720,7730) thus meaning 
that the full class of corresponding f{R) theories is viable. As a consequence, we may conclude that the present day 
SNela data are unable to discriminate between a A dominated Universe and this class of fourth order gravity theories. 

As a next step, we consider a SNAP -like survey [3| thus setting {AfsNeia, crm) — (2000,0.02). We use the same 
redshift distribution in Table 1 of and add 300 nearby SNela in the redshift range (0.03, 0.08). The Fisher matrix 
calculation gives for the uncertainties on the cosmographic parameters : 

(ai,(72,(T3,a4) - (0.08,1.0,4.8,13.7) . 

The significant improvement of the accuracy in the determination of {qo,jo,soJo) translates in a reduction of the 
errors on (7720, 7730) which now read : 

(^720,^^30) = (0.007,0.008) 

having assumed that, when SNAP data will be available, the matter density parameter Qm has been determined with 
a precision gm/^m ~ 1%- Looking again at Figs. [1] and O it is clear that the situation is improved. Indeed, the 
constraints on 7^20 makes it possible to narrow the range of allowed models with low matter content (the dashed line) , 
while models with typical values of VIm are still viable for wq covering almost the full horizontal axis. On the other 
hand, the constraint on 7730 is still too weak so that almost the full region plotted is allowed. 

Finally, we consider an hypothetical future SNela survey working at the same photometric accuracy as SNAP and 
with the same redshift distribution, but increasing the number of SNela up to AfsNeia — 6x10** as expected from, 
e.g., DES fl^, PanSTARRS SKYMAPPER ^4^, while still larger numbers may potentially be achieved by 
ALPACA [5ll| and LSST |l6l| . Such a survey can achieve : 

(CTi,(T2,CT3,(T4) = (0.02,0.2,0.9,2.7) 

so that, with (Jm/^m ~ 0.1%, we get : 

(^T20,(T3o) = (0.0015,0.0016) . 

Fig.[T]shows that, with such a precision on 7720, the region of wq values allowed essentially reduces to the ACDM value, 
while, from Fig. [21 it is clear that the constraint on 7730 definitively excludes models with low matter content further 
reducing the range of wq values to quite small deviations from the wq = —1. We can therefore conclude that such a 
survey will be able to discriminate between the concordance ACDM model and all the f{R) theories giving the same 
cosmographic parameters as quiessence models other than the ACDM itself. 

A similar discussion may be repeated for f{R) models sharing the same {qo,jo, so, ^o) values as the CPL model even 
if it is less intuitive to grasp the efficacy of the survey being the parameter space multivalued. For the same reason, 
we have not explored what is the accuracy on the double power - law or HS models, even if this is technically possible. 
Actually, one should first estimate the errors on the present day value of f{R) and its three time derivatives and then 
propagate them on the model parameters using the expressions obtained in Sect. VI. The multiparameter space to 
be explored makes this exercise quite cumbersome so that we leave it for a furthcoming work where we will explore 
in detail how these models compare to the present and future data. 
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VII. WHAT WE HAVE LEARNT FROM COSMOGRAPHY 

The recent amount of good quality data have given a new input to the observational cosmology. As often in science, 
new and better data lead to unexpected discoveries as in the case of the nowadays accepted evidence for cosmic 
acceleration. However, a fierce and strong debate is still open on what this cosmic speed up impHes for theoretical 
cosmology. The equally impressive amount of different (more or less) viable candidates have also generated a great 
confusion so that model independent analyses are welcome. A possible solution could come from cosmography rather 
than assuming ad hoc solutions of the cosmological Friedmann equations. Present day and future SNela surveys have 
renewed the interest in the determination of the cosmographic parameters so that it is worth investigating how these 
quantities can constrain cosmological models. 

Motivated by this consideration, in the framework of metric formulation of f{R) gravity, we have here derived the 
expressions of the present day values of /(i?) and its first three derivatives as function of the matter density parameter 
i^M, the Hubble constant Hq and the cosmographic parameters {qo,jo, sqJq). Although based on a third order Taylor 
expansion of f{R), we have shown that such relations hold for a quite large class of models so that they are valid 
tools to look for viable f{R) models without the need of solving the mathematically difficult nonlinear fourth order 
differential field equations. 

Notwithstanding the common claim that we live in the era of precision cosmology, the constraints on {qo,ja, so,lo) 
are still too weak to efficiently apply the program we have outlined above. As such, we have shown how it is possible 
to estabHsh a Hnk between the popular CPL parameterization of the dark energy equation of state and the derivatives 
of f{R), imposing that they share the same values of the cosmographic parameters. This analysis has lead to the 
quite interesting conclusion that the only f{R) function able to give the same values of {qo, jo, so,lo) as the ACDM 
model is indeed f{R) = R + 2A. If future observations will tell us that the cosmographic parameters are those of 
the ACDM model, we can therefore rule out all f{R) theories satisfying the hypotheses underlying our derivation of 
Eqs. (|87l) - 1(90|) . Actually, such a result should not be considered as a no way out for higher order gravity. Indeed, 
one could still work out a model with null values of /"(i?o) and /"'(i?o) as required by the above constraints, but 
non - vanishing higher order derivatives. One could well argue that such a contrived model could be rejected on the 
basis of the Occam razor, but nothing prevents from still taking it into account if it turns out to be both in agreement 
with the data and theoretically well founded. 

If new SNela surveys will determine the cosmographic parameters with good accuracy, acceptable constraints on 
the two dimensionless ratios 7720 oc f"{Ro)/ f{Ro) and 7730 cx f"'{Ro)/f{Ro) could be obtained thus allowing to 
discriminate among rival f{R) theories. To investigate whether such a program is feasible, we have pursued a Fisher 
matrix based forecasts of the accuracy future SNela surveys can achieve on the cosmographic parameters and hence 
on (?720i?73o)- It turns out that a SNAP -like survey can start giving interesting (yet still weak) constraints allowing 
to reject f{R) models with low matter content, while a definitive improvement is achievable with future SNela survey 
observing ~ 10'* objects thus making it possible to discriminate between ACDM and a large class of fourth order 
theories. It is worth stressing, however, that the measurement of should come out as the result of a model 
independent probe such as the gas mass fraction in galaxy clusters which, at present, is still far from the 1% requested 
precision. On the other hand, one can also rely on the flM estimate from the CMBR anisotropy and polarization 
spectra even if this comes to the price of assuming that the physics at recombination is strictly described by GR 
so that one has to limit its attention to f{R) models reducing to f{R) cc R during that epoch. However, such an 
assumption is quite common in many f{R) models available in literature so that it is not a too restrictive limitation. 

A further remark is in order concerning what kind of data can be used to constrain the cosmographic parameters. 
The use of the fifth order Taylor expansion of the scale factor makes it possible to not specify any underlying physical 
model thus relying on the minimalist assumption that the Universe is described by the flat Robertson - Walker metric. 
While useful from a theoretical perspective, such a generality puts severe limitations to the dataset one can use. 
Actually, we can only resort to observational tests depending only on the background evolution so that the range of 
astrophysical probes reduces to standard candles (such as SNela and possibly Gamma Ray Bursts |34]) and standard 
rods (such as the angular size - redshift relation for compact radiosources). Moreover, pushing the Hubble diagram 
to z ~ 2 may rise the question of the impact of gravitational lensing amplification on the apparent magnitude of the 
adopted standard candle. The magnification probability distribution function depends on the growth of perturbations 
[49I. I75I. [s^ . [stI . [sol so that one should worry about the underlying physical model in order to es tima t e whethe r this 
effect biases the estimate of the cosmographic parameters. However, it has been shown [sil . Isi . IIITI . Il34l . Il39l | that 
the gravitational lensing amplification does not alter significantly the measured distance modulus for z ~ 1 SNela. 
Although such an analysis has been done for GR based models, we can argue that, whatever is the f{R) model, the 
growth of perturbations finally leads to a distribution of structures along the line of sight that is as similar as possible 
to the observed one so that the lensing amplification is approximately the same. We can therefore argue that the 
systematic error made by neglecting lensing magnification is lower than the statistical ones expected by the future 
SNela surveys. On the other hand, one can also try further reducing this possible bias using the method of fiux 
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averaging |168l | even if, in such a case, our Fisher matrix calculation should be repeated accordingly. It is also worth 
noting that the constraints on the cosmographic parameters may be tightened by imposing some physically motivated 
priors in the parameter space. For instance, we can impose that the Hubble parameter H{z) stays always positive 
over the full range probed by the data or that the transition from past deceleration to present acceleration takes 
place over the range probed by the data (so that we can detect it). Such priors should be included in the likelihood 
definition so that the Fisher matrix should be recomputed which is left for a forthcoming work. 

Although the present day data are still too Hmited to efficiently discriminate among rival f{R) models, we are 
confident that an aggressive strategy aiming at a very precise determination of the cosmographic parameters could 
offer stringent constraints on higher order gravity without the need of solving the field equations or addressing the 
complicated problems related to the growth of perturbations. Almost 80 years after the pioneering distance - redshift 
diagram by Hubble, the old cosmographic approach appears nowadays as a precious observational tool to investigate 
the new developments of cosmology. 



VIII. THE WEAK-FIELD LIMIT OF /(i?)-GRAVITY 

Before facing the problem of galaxy clusters by /(i?)-gravity, a discussion is due on the weak-field limit of such a 
theory which, being of fourth order in metric formalism, could lead to results radically different with respect to the 
case f{R) = R, the standard second order General Relativity. 

Let us consider the general action : 



A - 



J d^x^[f(R) + X£„,] 



(124) 



WttG 



is the 



where f{R) is an analytic function of the Ricci scalar R, g is the determinant of the metric g^j^, X 

coupling constant and Cm is the standard perfect-fluid matter Lagrangian. Such an action is the straightforward 
generalization of the Hilbert-Einstein action of GR obtained for f{R) — R. Since we are considering the metric 
approach, field equations are obtained by varying l|124p with respect to the metric : 



1 A! 



(125) 



where T,, 



_ -2 Sj^C^) 



^i^f — — — the energy momentum tensor of matter, the prime indicates the derivative with respect 
to R and □ — . We adopt the signature (+,—,—,—). 

As discussed in details in [sO], we deal with the Newtonian and the post-Newtonian Hmit of f{R) - gravity on 
a spherically symmetric background. Solutions for the field equations can be obtained by imposing the spherical 
symmetry [29l |: 



goo (x" , r)dx^ ^ -I- g^r {x^ , r)dr'^ — d^l 



(126) 



where — ct and dfl is the angular element. 

To develop the post-Newtonian limit of the theory, one can consider a perturbed metric with respect to a Minkowski 
background .g^,y — 77^,^ -t- hfj^^. The metric coefficients can be developed as: 



f gu{t,r)c^l+gif{t,r)+git\t,r) 
grr{t,r) ~ -1 + gil\t,r) 



(127) 



-r^ sin^ I 



where we put, for the sake of simplicity, c — 1 , x° — ct ^ t. We want to obtain the most general result without 
imposing particular forms for the /(i?)-Lagrangian. We only consider analytic Taylor expandable functions 



/(i?)~/0 + /li?,+ /2i?'+/3i?' 



(128) 
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To obtain the post-Newtonian approximation of f{R) - gravity, one has to plug the expansions l|127p and (|128p into 
the field equations (|125p and then expand the system up to the orders 0(0), 0(2) and 0(4) . This approach provides 
general results and specific (analytic) Lagrangians are selected by the coefficients fi in l|128p [30]. 

If we now consider the 0(2) - order of approximation, the field equations l|125p . in the vacuum case, results to be 



(2) 



f (2) 
' Jl''^9tt,rr 



= 



-fi9tti - fi9rl]r + W^? + 4/2ri?: 
/iri?(2) + 6/2[2i?,l'^ + rR^pr] = 



(2)i 



(129) 



2(?^^^ + r[2g^u]r " ^^^^^ + '^9?r,r 



(2) 1 







It is evident that the trace equation (the fourth in the system (|129p ). provides a differential equation with respect to 
the Ricci scalar which allows to solve the system at 0(2) - order. One obtains the general solution : 
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e^-'r 



h 
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R. 



where ^ = — -, /i and /2 are the expansion coefficients obtained by the /(i?)-Taylor series. In the Hmit / 

6/2 

for a point-Hke source of mass M we recover the standard Schwarzschild solution. Let us notice that the integration 
constant 5q is dimensionless, while the two arbitrary time-functions 5i{t) and 52{t) have respectively the dimensions of 
lenght~^ and lenght~^; ^ has the dimension lenght~'^. As extensively discussed in [sOl, the functions 6i{t) (i = 1,2) 
are completely arbitrary since the differential equation system (|129p depends only on spatial derivatives. Besides, the 
integration constant 6o can be set to zero, as in the standard theory of potential, since it represents an unessential 
additive quantity. In order to obtain the physical prescription of the asymptotic fiatness at infinity, we can discard 
the Yukawa growing mode in (|130p and then the metric is : 



The Ricci scalar curvature is 



2GA/ Si{t)e- 
2GA/ 



fir 
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The solution can be given also in terms of gravitational potential. In particular, we have an explicit Newtonian-like 
term into the definition. The first of (|13Qp provides the second order solution in term of the metric expansion (see the 
definition (|127p ). In particular, it is gtt = 1 ' — t i ^'^ > 

/(i?)-theory is 



'24)grav = 1 -f git and then the gravitational potential of an analytic 



GM Si{t)e- 



fir 



6^r 



(133) 
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Among the possible analytic /(i?)-models, let us consider the Taylor expansion where the cosmological term (the 
above /o) and terms higher than second have been discarded. For the sake of simplicity, we rewrite the Lagrangian 
(fT28l) as 

J{R)^ aiR + a2R^ + - (134) 

and specify the above gravitational potential (|133p . generated by a point-like matter distribution, as: 

iGM / 1 , , ^ 

1 + ^e-T: , 135 



' ' Aair V 3 
where 

L = L(ai,a2)- f-^j . (136) 

L can be defined as the interaction length of the problem^^ due to the correction to the Newtonian potential. We 
have changed the notation to remark that we are doing only a specific choice in the wide class of potentials (|133p , but 
the following considerations are completely general. 

IX. EXTENDED SYSTEMS 

The gravitational potential (|135p is a point-like one. Now we have to generalize this solution for extended systems. 
Let us describe galaxy clusters as spherically symmetric systems and then we have to extend the above considerations 
to this geometrical configuration. We simply consider the system composed by many infinitesimal mass elements dm 
each one contributing with a point-like gravitational potential. Then, summing up all terms, namely integrating them 
on a spherical volume, we obtain a suitable potential. Specifically, we have to solve the integral: 

<^{r)^ r'^dr' i sinO' d9' duu' 4>{r') . (137) 
Jo ^0 Ja 

The point-like potential (|135p can be split in two terms. The Newtonian component is 

3GM , , 

Mr)=-^ (138) 

The extended integral of such a part is the well-known (apart from the numerical constant ^) expression. It is 

^ , , 3 GM(< r) 

$M.) = -^^^ (139) 

where M (< r) is the mass enclosed in a sphere with radius r. The correction term: 

, , , GM e-i 

Mr)-~^— (140) 
considering some analytical steps in the integration of the angular part, gives the expression: 

2, — I n QQ _ I ' ^1 _ I ^ I 

<^c(r) L / dr'r'pir') 141 

The radial integral is numerically estimated once the mass density is given. We underline a fundamental difference 
between such a term and the Newtonian one: while in the latter, the matter outside the spherical shell of radius r 
does not contribute to the potential, in the former external matter takes part to the integration procedure. For this 
reason we split the corrective potential in two terms: 



Such a length is function of the series coefficients, ai and 02, and it is not a free independent parameter in the following fit procedure. 
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if r' < r: 



r _ \-r-r'\ 

^C,^■at(r)^ J- ■ L / dr'r'p{r') 

L / dr r p{r ) ■ e 



4 







if r' > 



r. ^ ,.00 '■+■■'1 

2nG ^ f , , , , e — - e — 



^cMr) = -^-Lj dr'r'pir') 



r 



• L / dr r p{r ) ■ e ^ 



The total potential of the spherical mass distribution will be 

$(r) = $Ar(r) + $a»«t(r-) + $c,ext(r) (142) 

As we will show below, for our purpose, we need the gravitational potential derivative with respect to the variable 
r; the two derivatives may not be evaluated analytically so we estimate them numerically, once we have given an 
expression for the total mass density p{r). While the Newtonian term gives the simple expression: 

- = ^ (143) 

dr 4ai 

The internal and external derivatives of the corrective potential terms are much longer. We do not give them explicitly 
for sake of brevity, but they are integral-functions of the form 

fP(r) 

T{r,r') = / dr' f{r,r) (144) 

J a(r) 



from which one has: 

dr dr 

-/(r,a(r))^(r) + /(r,/3(r))^(r) (145) 

Such an expression is numerically derived once the integration extremes are given. A general consideration is in order 
at this point. Clearly, the Gauss theorem holds only for the Newtonian part since, for this term, the force law scales 
as 1/r^. For the total potential (|135p . it does not hold anymore due to the correction. From a physical point of view, 
this is not a problem because the full conservation laws are determined, for /(i^) -gravity, by the contracted Bianchi 
identities which assure the self-consistency. For a detailed discussion, see [IB, [Z2| ■ 



X. THE CLUSTER MASS PROFILES 



Clusters of galaxies are generally considered self-bound gravitational systems with spherical symmetry and in 
hydrostatic equilibrium if virialized. The last two hypothesis are still widely used, despite of the fact that it has been 
widely proved that most clusters show more complex morphologies and/or signs of strong interactions or dynamical 
activity, especially in their innermost regions (|42.[5a|'). 

Under the hypothesis of spherical symmetry in hydrostatic equihbrium, the structure equation can be derived from 
the collisionless Boltzmann equation 

UpoUr) al) + ^-P^ia'^ „ = -p.^r) ■ ^ (146) 
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where $ is the gravitational potential of the cluster, cr^ and erg are the mass-weighted velocity dispersions in the 
radial and tangential directions, respectively, and p is gas mass-density. For an isotropic system, it is 



The pressure profile can be related to these quantities by 

P(r) = alpgasir) 

Substituting Eqs. (|147p and l|148p into Eq. (|146p . we have, for an isotropic sphere. 



dP{r) 



-Pgas{r) 



dr ' dr 

For a gas sphere with temperature profile T(r), the velocity dispersion becomes 

2 kT{r) 

(7,. = 



(147) 
(148) 

(149) 
(150) 



where k is the Boltzmann constant, fi ~ 0.609 is the mean mass particle and nip is the proton mass. Substituting 
Eqs. lfT48l) and ifTSOj) into Eq. (fT49| . we obtain 



d (kT{r) \ d$ 



or, equivalently. 



d$ kT{r) 



dlnr 



dr fiiripr 

Now the total gravitational potential of the cluster is: 

$(r) = $Af(r) -I- $c(r) 

with 



c?lnr 



(151) 



(152) 



(153) 



It is worth underlining that if we consider only the standard Newtonian potential, the total cluster mass Mci,N{r) is 
composed by gas mass -I- mass of galaxies -I- cD-galaxy mass -I- dark matter and it is given by the expression: 



Mcl,N{r)^Mgasir) + Mgal{r) + McDgal{r) + MDM(r) 



kT{r) 
p,mpG 



d\npgas{r) dlnT{r) 



dlnr 



dlnr 



(154) 



Mci,N means the standard estimated Newtonian mass. Generally the galaxy part contribution is considered negligible 
with respect to the other two components so we have: 



Mci,N{r) « Mgas{r) + MoAiir) 



kT(r) 

1 

pirip 



dlnpgas{r) d\nT{r) 



dlni 



In) 



Since the gas-mass estimates are provided by X-ray observations, the equilibrium equation can be used to derive the 
amount of dark matter present in a cluster of galaxies and its spatial distribution. 

Inserting the previously defined extended- corrected potential of Eq. (|152p into Eq. (|15ip . we obtain: 
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dr 
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(155) 
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from which the extended- corrected mass estimate follows: 



/ N 4ai , d^c" / s 



4ai 
~3~ 



fcT(r) / dlnpgasir) d\nT{r) 

r \ 

^nipG \ dlnr dlnr 



(156) 



Since the use of a corrected potential avoids, in principle, the additional requirement of dark matter, the total cluster 
mass, in this case, is given by: 



Mci,Ec{r) = Mgasir) + Mgal{r) + McDgal{r) 

and the mass density in the <i>c term is 

Pcl,Ec{r) = Pgasir) + pgal{r) + PCDgal{r) 



(157) 



(158) 



with the density components derived from observations. 

In this work, we will use Eq. I|156p to compare the baryonic mass profile Md^Ecif), estimated from observations, 

with the theoretical deviation from the Newtonian gravitational potential, given by the expression — — ^-^"^ , {^)- 
Our goal is to reproduce the observed mass profiles for a sample of galaxy clusters. 



XI. THE GALAXY CLUSTER SAMPLE 



The form alism des cribed in § [X] can be applied to a sample of 12 galaxy clusters. We shall use the cluster sample 
studied in [l64L Il65l | which consists of 13 low-redshift clusters spanning a temperature range 0.7 9.0 keV derived 
from high quality Chandra archival data. In all these clusters, the surface brightness and the gas temperature profiles 
are measured out to large radii, so that mass estimates can be extended up to r5oo or beyond. 



A. The Gas Density Model 

The gas density distribution of the clusters in the sample is described by the analytic model proposed in |165 
Such a model modifies the classical /3— model to represent the characteristic properties of the observed X-ray surface 
brightness profiles, i.e. the power-law-type cusps of gas density in the cluster center, instead of a fiat core and the 
steepening of the brightness profiles at large radii. Eventually, a second /3— model, with a small core radius, is added 
to improve the model close to the cluster cores. The analytical form for the particle emission is given by: 



(r/re)-" 1 



(l + r2/7-2)3'3-"/2 [l + r^/riyh 



[l + r-^/rl^fP^ 

which can be easily converted to a mass density using the relation: 



(159) 



1.4 , , 

Pgas = riT ■ finip = —Hemp (160) 

where nx is the total number density of particles in the gas. The resulting model has a large number of parameters, 
some of which do not have a direct physical interpretation. While this can often be inappropriate and computationally 
inconvenient, it suits well our case, where the main requirement is a detailed qualitative description of the cluster 
pro files . 

In |l65l |. Eq. I|159p is applied to a restricted range of distances from the cluster center, i.e. between an inner cutoff 
rmin, chosen to exclude the central temperature bin (« 10 20 kpc) where the ICM is Hkely to be multi-phase, and 
rdet, where the X-ray surface brightness is at least 3a significant. We have extrapolated the above function to values 
outside this restricted range using the following criteria: 

• for r < rmin, we have performed a linear extrapolation of the first three terms out to r = kpc; 
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Table I: Column 1: Cluster name. Column2: Richness. Column 3: cluster total mass. Column 4: gas mass. Column 5: galaxy 
mass. Column 6: cD-galaxy mass. All mass values are estimated at r = rmax- Column 7: ratio of total galaxy mass to gas 
mass. Column 8: minimum radius. Column 9: maximum radius. 



name 


K 




Mgas 


Mgal 


McDgal 


gal 
gas 










(Mo) 


(Mo 


) 


(Mo 


) 


(Mo) 




(kpc) (kpc) 


A133 





4.35874 ■ 


10" 


2.73866 • 


10^3 


5.20269 • 


10^2 


1.10568 • 


10^2 


0.23 


86 


1060 


A262 





4.45081 ■ 


10^3 


2.76659 ■ 


10^2 


1.71305 • 


10" 


5.16382 • 


10^2 


0.25 


61 


316 


A383 


2 


2.79785 ■ 


10" 


2.82467 ■ 


10^3 


5.88048 • 


10^2 


1.09217 • 


10^2 


0.25 


52 


751 


A478 


2 


8.51832 ■ 


10" 


1.05583 • 


10" 


2.15567 • 


10^3 


1.67513 • 


10^2 


0.22 


59 


1580 


A907 


1 


4.87657 ■ 


10" 


6.38070 • 


10^3 


1.34129 • 


10^3 


1.66533 ■ 


10^^ 


0.24 


563 


1226 


A1413 


3 


1.09598 ■ 


10^5 


9.32466 • 


10^3 


2.30728 • 


10^3 


1.67345 • 


10^2 


0.26 


57 


1506 


A1795 


2 


5.44761 ■ 


10" 


5.56245 ■ 


10^3 


4.23211 • 


10^2 


1.93957 • 


10^2 


0.11 


79 


1151 


A1991 


1 


1.24313 ■ 


10" 


1.00530 • 


10^3 


1.24608 • 


10^^ 


1.08241 ■ 


10^2 


0.23 


55 


618 


A2029 


2 


8.92392 ■ 


10" 


1.24129 • 


10" 


3.21543 • 


10^3 


1.11921 • 


10^^ 


0.27 


62 


1771 


A2390 


1 


2.09710 ■ 


10^^ 


2.15726 • 


10" 


4.91580 • 


10^3 


1.12141 ■ 


10^^ 


0.23 


83 


1984 


MKW4 




4.69503 ■ 


10^3 


2.83207 • 


10^2 


1.71153 • 


10" 


5.29855 ■ 


10" 


0.25 


60 


434 


RXJ1159 




8.97997 ■ 


10^3 


4.33256 ■ 


10^2 


7.34414 • 


10" 


5.38799 ■ 


10" 


0.29 


64 


568 



• for r > rdet, we have performed a linear extrapolation of the last three terms out to a distance f for which 
Pgasir) = Pc Pc being the critical density of the Universe at the cluster redshift: pc — Pc,Q ■ (1 + z)^- For radii 
larger than f, the gas density is assumed constant at Pgas(r). 

We point out t hat, i n Table [H the radius limit Tmin is almost the same as given in the previous definition. When the 
value given by |l65l | is less than the cD-galaxy radius, which is defined in the next section, we choose this last one 
as the lower limit. On the contrary, Vmax is quite different from Vdet- it is fixed by considering the higher value of 
temperature profile and not by imaging methods. 

We then compute the gas mass Mgas{r) and the total mass Mci,N{r), respectively, for all clusters in our sample, 
substituting Eq. I|159p into Eqs. i|160p and l|154p . respectively; the gas temperature profile has been described in 
details in § IXIBI The resulting mass values, estimated at r = r^ax, are listed in Table [H 



B. The Temperature Profiles 

As stressed in S IXI A[ for the purpose of this work, we need an accurate qualitative description of the radial beha vior 
of the gas properties. Standard isothermal or polytropic models, or even the more complex one proposed in [l65l |. do 
not provide a good description of the data at all radii and for all clusters in the present sample. We hence describe 
the gas temperature profiles using the straightforward X-ray spectral analysis results, without the introduction of any 
analytic model. 

X-ray spectral values have been provided by A. Vikhlinin (private communication). A detailed description of the 
relative spectral analysis can be found in [l6J|. 



C. The Galaxy Distribution Model 



The galaxy density can be modelled as proposed by Even if the galaxy distribution is a pomi-distribution 
instead of a continuous function, assuming that galaxies are in equilibrium with gas, we can use a model, oc r~^, 
for r < Rc from the cluster center, and a steeper one, oc r~^-^, for r > Rc, where Rc is the cluster core radius (its 
value is taken from Vikhlinin 2006). Its final expression is: 



Pgaiir) = 



qal.l 



ga/,2 • 



r <Rr. 



> Rc. 



(161) 
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Table II: Column 1: Cluster name. Column 2: first derivative coefRcient, ai, of f(R) series. ColumnS: Icr confidence interval for ai. 
Column 4: second derivative coefficient, a2, of f(R) series. Column 5: Itr confidence interval for 02. Column 6: characteristic length, L, 
of the modified gravitational potential, derived from ai and 02- Column 7 : Icr confidence interval for L. 



name 




ai 


[ai - 


- If^, 


ai + Icr] 


£I2 


]a2 — lo", a-i + Icr] 


L 


[L — Icr, L + Icr] 














(kpc ) 


/I 2\ 

(kpc ) 


(kpc) 


(kpc) 


A133 


0, 


,085 


[0, 


,078, 


0.091] 


-4.98 ■ 10^ 


]-2.38 ■ 10^, -1.38 ■ 10^] 


591.78 


]323.34, 1259.50] 


A262 


0, 


,065 


[0, 


,061, 


0.071] 


-10.63 


]-57.65, -3.17] 


31.40 


[17.28, 71.10] 


A383 


0, 


,099 


[0, 


,093, 


0.108] 


-9.01 ■ 10^ 


]-4.10 ■ 10^ -3.14 • 10^] 


234.13 


]142.10, 478.06] 


A478 


0, 


,117 


[0, 


,114, 


0.122] 


-4.61 ■ 10^ 


]-1.01 ■ 10"*, -2.51 • 10^] 


484.83 


]363.29, 707.73] 


A907 


0, 


,129 


[0, 


,125, 


0.136] 


-5.77- 10^ 


]-1.54 ■ 10"*, -2.83 ■ 10^] 


517.30 


]368.84, 825.00] 


A1413 


0, 


,115 


[0, 


,110, 


0.119] 


-9.45 ■ lO"" 


]-4.26 ■ 10^ -3.46 • 10"] 


2224.57 


[1365.40, 4681.21] 


A1795 


0, 


,093 


[0, 


,084, 


0.103] 


-1.54- 10^ 


]-1.01 • lO'*, -2.49 ■ 10^] 


315.44 


]133.31, 769.17] 


A1991 


0, 


,074 


[0, 


,072, 


0.081] 


-50.69 


]-3.42 ■ 10^ -13] 


64.00 


132.63, 159.40] 


A2029 


0, 


,129 


[0, 


,123, 


0.134] 


-2.10 ■ 10" 


]-7.95-10'*, -8.44-10^] 


988.85 


]637.71, 1890.07] 


A2390 


0, 


,149 


[0, 


,146, 


0.152] 


-1.40 ■ 10" 


]-5.71•10^ -4.46-10'^] 


7490.80 ]4245.74, 15715.60] 


MKW4 


0, 


,054 


[0, 


,049, 


0.060] 


-23.63 


]-1.15 ■ 10^ -8.13] 


51.31 


]30.44, 110.68] 


RXJ1159 


0, 


,048 


[0, 


,047, 


0.052] 


-18.33 


]-1.35 ■ 10^ -4.18] 


47.72 


]22.86, 125.96] 



where the constants Pgai,i and Pgai,2 are chosen in the following way: 

• [il provides the central number density of galaxies in rich compact clusters for galaxies located within a 1.5 
h~^Mpc radius from the cluster center and brighter than + 2™ (where is the magnitude of the third 
brightest galaxy): nga/,o ~ lO^ft'^ galaxies Mpc""^. Then we fix Pgai,i in the range ^ 10'^'* 10^^ kg/kpc^. For 
any cluster obeying the condition chosen for the mass ratio gal-to-gas, we assume a typical elliptical and cD 
galaxy mass in the range 10^^ 10^^ Mq. 

• the constant Pgai,2 has been fixed with the only requirement that the galaxy density function has to be continuous 
at Re- 

We have tested the effect of varying galaxy density in the above range ^ 10'^^ 10"^® kg/kpc'^ on the cluster with the 
lowest mass, namely A262. In this case, we would expect great variations with respect to other clusters; the result is 
that the contribution due to galaxies and cD-galaxy gives a variation < 1% to the final estimate of fit parameters. 
The cD galaxy density has been modelled as described in [l43| ; they use a Jaffe model of the form: 



PCDgal — 




where rc is the core radius while the central density is obtained from Mj = -TrR^po,j. The mass of the cD galaxy 

has been fixed at 1.14 x 10^^ Mq, with rc = Re/0.76, with Re = 25 kpc being the effective radius of the galaxy. The 
central galaxy for each cluster in the sample is assumed to have approximately this stellar mass. 

We have assumed that the total galaxy-component mass (galaxies plus cD galaxy masses) is ~ 20 25% of the 
gas mass: in |l4Cll |. the mean fraction of gas versus the total mass (with dark matter) for a cluster is estimated to be 
15 ~ 20%, while the same quantity for galaxies is 3 5%. This means that the relative mean mass ratio gal-to-gas 
in a cluster is ~ 20 25%. We have varied the parameters Pgai,i, Pgai,2 and Mj in their previous defined ranges to 
obtain a mass ratio between total galaxy mass and total gas mass which lies in this range. Resulting galaxy mass 

values and ratios — — , estimated at r = rmax , are listed in Table [H 

gas 

In Fig. (1), we show how each component is spatially distributed. The CD-galaxy is dominant with respect to the 
other galaxies only in the inner region (below 100 kpc). As already stated in S IXI A[ cluster innermost regions have 
been excluded from our analysis and so the contribution due to the cD-galaxy is practically negligible in our analysis. 
The gas is, as a consequence, clearly the dominant visible component, starting from innermost regions out to large 
radii, being galaxy mass only 20 25% of gas mass. A similar behavior is shown by all the clusters considered in our 
sample. 
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D. Uncertainties on mass profiles 

Uncertainties on the cluster total mass profiles have been estimated performing Monte-Carlo simulations |l08l |. We 
proceed to simulate temperature profiles a nd c hoose random radius-temperature values couples for each bin which 
we have in our temperature data given by |l64| . Random temperature values have been extracted from a Gaussian 
distribution centered on the spectral values, and with a dispersion fixed to its 68% confidence level. For the radius, 
we choose a random value inside each bin. We have performed 2000 simulations for each cluster and perform two 
cuts on the simulated profile. First, we exclude those profiles that give an unphysical negative estimate of the mass: 
this is possible when our simulated couples of quantities give rise to too high temperature-gradient. After this cut, 
we have « 1500 simulations for any cluster. Then we have ordered the resulting mass values for increasing radius 
values. Extreme mass estimates (outside the 10 90% range) are excluded from the obtained distribution, in order 
to avoid other high mass gradients which give rise to masses too different from real data. The resulting limits provide 
the errors on the total mass. Uncertainties on the electron-density profiles has not been included in the simulations, 
being them negligible with respect to those of the gas-temperature profiles. 



E. Fitting the mass profiles 

In the above sections, we have shown that, with the aid of X-ray observations, modelling theoretically the galaxy 
distribution and using Eq. (|156p . we obtain an estimate of the baryonic content of clusters. 
We have hence performed a best-fit analysis of the theoretical Eq. (|156p 



Mbar,th{r)-- 



versus the observed mass contributions 



4a 1 
4ai 



kT{r) ( din Pgasjr) d In T(r) 



finipG 

d$c, 
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3G dr 



■(r) 
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Mbar,obsir) = Mgas{r) + Mgal{r) + McDgal{r) 



(164) 



Since not all the data involved in the above estimate have measurable errors, we cannot perform an exact x-square 
minimization: Actually, we can minimize the quantity: 



1 



N 



,obs 



Mbar,theoY 



N -n„ 



1 ^ 

1=1 



Mbar. 



the 



(165) 



where N is the number of data and rip = 2 the free parameters of the model. We minimize the x-square using 
the Markov Chain Monte Carlo Method (MCMC). For each cluster, we have run various chains to set the best 
parameters of the used algorithm, the Metropolis-Hastings one: starting from an initial parameter vector p (in our 
case p = (ai,02)), we generate a new trial point p' from a tested proposal density g(p',p), which represents the 
conditional probability to get p', given p. This new point is accepted with probability 



a(P,p') = 



£(d|pOP(pOg(p',p) 
L(d|p)P(p)g(p,p') 



where d are the data, L(d|p') oc exp(— x^/2) is the likelihood function, -P(p) is the prior on the parameters. In 
our case, the prior on the fit parameters is related to Eq. 



(|136p : being L a length, we need to force the ratio 
ai/a2 to be positive. The proposal density is Gaussian symmetric with respect of the two vectors p and p', namely 
g(p, p') cx exp(— Ap^/2cr^), with Ap = p — p'; we decide to fix the dispersion a of any trial distribution of parameters 
equal to 20% of trial oi and a2 at any step. This means that the parameter a reduces to the ratio between the 
likelihood functions. 

We have run one chain of 10^ points for every cluster; the convergence of the chains has been tested using the 
power spectrum analysis from [63]. The key idea of this method is, at the same time, simple and powerful: if we 
take the power spectra of the MCMC samples, we will have a great correlation on small scales but, when the chain 
reaches convergence, the spectrum becomes fiat (Hke a white noise spectrum); so that, by checking the spectrum of 
just one chain (instead of many parallel chains as in Gelmann-Rubin test) will be sufficient to assess the reached 
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Figure 5: Left panel: histogram of the sample points for parameter ai in Abell 383 coming out the MCMC implementation used to 
estimate best fit values and errors for our fitting procedure as described in § IXI El Binning (horizontal axis) and relative frequencies 
(vertical axis) are given by automatic procedure from Mathematica6.0. Right panel: power spectrum test on sample chain for parameter 
ai using the method described in S IXI El Black line is the logarithm of the analytical template Eq. 1)16811 for power spectrum; gray line is 
the discrete power spectrum obtained using Eq. 1II66II - II167II . 



convergence. Remanding to [61 
power spectrum of the chains: 
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for a detailed discussion of all the mathematical steps. Here we calculate the discrete 
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where N and Xn are the length and the element of the sample from the MCMC, respectively, j = 1, . . . — 1. The 
wavenumber kj of the spectrum is related to the index j by the relation kj = ^ . Then we fit it with the analytical 
template: 



P(fc) = Po 



or in the equivalent logarithmic form: 



In R = In Po + In 



{k* ikY 
1 + (fc*/fc)" 



\^(k*ik,r 



l + rj 



(168) 



(169) 



where 7 — 0.57216 is the Euler-Mascheroni number and rj are random measurement errors with < rj >— and 
< TiVj >= Sijn^/6. From the fit, we estimate the two fundamental parameters, Po and j* (the index corresponding 
to k*). The first one is the value of the power spectrum extrapolated for fc ^ and, from it, we can derive the 

convergence ratio from r « — ; if r < 0.01, we can assume that the convergence is reached. The second parameter is 

related to the turning point from a power-law to a fiat spectrum. It has to be > 20 in order to be sure that the number 
of points in the sample, coming from the convergence region, are more than the noise points. If these two conditions 
are verified for all the parameters, then the chain has reached the convergence and the statistics derived from MCMC 
well describes the underlying probability distribution (typical results are shown in Figs. (2)-(3)). Following [6^ 
prescriptions, we perform the fit over the range 1 < j < jmax, with jmax ^ lOj*, where a first estimation of j* 
can be obtained from a fit with jmax — 1000, and then performing a second iteration in order to have a better 
estimation of it. Even if the convergence is achieved after few thousand steps of the chain, we have decided to run 
longer chains of 10^ points to reduce the noise from the histograms and avoid under- or over- estimations of errors 
on the parameters. The i — a confidence levels are easily estimated deriving them from the final sample the 15.87-th 
and 84.13-th quantiles (which define the 68% confidence interval) for i = 1, the 2.28-th and 97. 72-th quantiles (which 
define the 95% confidence interval) for z = 2 and the 0.13-th and 99.87-th quantiles (which define the 99% confidence 
interval) for z = 3. 

After the description of the method, let us now comment on the achieved results. 
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XII. RESULTS 

The numerical results of our fitting analysis are summarized in Table 2; we give the best fit values of the independent 
fitting parameters ai and 02, and of the gravitational length L, considered as a function of the previous two quantities. 
In Figs. (3)- (5), we give the typical results of fitting, with histograms and power spectrum of samples derived by the 
MCMC, to assess the reached convergence (fiat spectrum at large scales). 

The goodness and the properties of the fits are shown in Figs. (6)- (17). The main property of our results is the 
presence of a typical scale for each cluster above which our model works really good (typical relative differences are 
less than 5%), while for lower scale there is a great difference. It is possible to see, by a rapid inspection, that this 
turning-point is located at a radius « 150 kpc. Except for very large clusters, it is clear that this value is independent 
of the cluster, being approximately the same for any member of the considered sample. 

There are two main independent explanations that could justify this trend: Hmits due to a break in the state of 
hydrostatic equilibrium or limits in the series expansion of the /(i?)-models. 

If the hypothesis of hydrostatic equilibrium is not co rrect , then we are in a regime where the fundamental relations 
Eqs. (|146p - (|15ip . are not working. As discussed in [l64l |. the central (70 kpc) region of every cluster is strongly 
affected by radiative cooling and thus it cannot directly be related to the depth of the cluster potential well. This 
means that, in this region, the gas is not in hydrostatic equilibrium but in a multi-phase, turbulent state, mainly 
driven by some astrophysical, non-gravitational interaction. In this case, the gas cannot be used as a good standard 
tracer. 

We have also to consider another limit of our modeUing: the requirement that the /(i?)-function is Taylor ex- 
pandable. The corrected gravitational potential which we have considered is derived in the weak field limit, which 
means 

R-Ro«— (170) 

where Rq is the background value of the curvature. If this condition is not satisfied, the approach does not work 
(see js^l for a detailed discussion of this point). Considering that ai/a2 has the dimension of length~'^ this condition 
defines the length scale where our series approximation can work. In other words, this indicates the Hmit in which 
the model can be compared with data. 

For the considered sample, the fit of the parameters ai and 02, spans the length range {19; 200} kpc (except for the 
biggest cluster). It is evident that every galaxy cluster has a proper gravitational length scale. It is worth noticing 
that a similar situation, but at completely different scales, has been found out for low surface brightness galaxies 
modelled by /(i?)-gravity (2^ . 

Considering the data at our disposal and the analysis which we have performed, it is not possible to quantify exactly 
the quantitative amount of these two different phenomena (i.e. the radiative cooling and the validity of the weak 
field Hmit). However, they are not mutually exclusive but should be considered in details in view of a more refined 
modelling 

Similar issues are present also in [l^: they use the the Metric - Skew - Tensor - Gravity (MSTG) as a generalization 
of the Einstein General Relativity and derive the gas mass profile of a sample of clusters with gas being the only 
baryonic component of the clusters. They consider some clusters included in our sample (in particular, A133, A262, 
A478, A1413, A1795, A2029, MKW4) and they find the same different trend for r < 200 kpc, even if with a different 
behavior with respect of us: our model gives lower values than X-ray gas mass data while their model gives higher 
values with respect to X-ray gas mass data. This stresses the need for a more accurate modelling of the gravitational 
potential. 

However, our goal is to show that potential (|135p is suitable to fit the mass profile of galaxy clusters and that it 
comes from a self-consistent theory. 

In gen eral, it can be shown that the weak field limit of extended theories of gravity has Yukawa-like corrections 
[gil . Il52ll. Specificallv. given theory of gravity of order {2n + 2), the Yukawa corrections to the Newtonian potential 
are n |129l |. This means that if the effective Lagrangian of the theory is 

c = f{R,aR,.n''R,.n'^R)y/^ (171) 



^ other secondary phenomena as cooling flows, merger and asymmetric shapes have to be considered in view of a detailed modelling of 
clusters. However, in this work, we are only interested to show that extended gravity could be a valid alternative to dark matter in 
order to explain the cluster dynamics. 
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Standard General Relativity, where Yukawa corrections are not present, is recovered for n = (second order theory) 
while the /(i?)-gravity is obtained for n = 1 (fourth-order theory). Any □ operator introduces two further derivation 
orders in the field equations. This kind of Lagrangian comes out when quantum field theory is formulated on curved 
spacetime In the series (|172p . G is the value of the gravitational constant considered at infinity, L/c is the 

interaction length of the k-th component of the non-Newtonian corrections. The amplitude ak of each component is 
norm alized to the standard Newtonian term; the sign of ak tells us if the corrections are attractive or repulsive (see 
[l7ll | for details). Moreover, the variation of the gravitational coupHng is involved. In our case, we are taking into 
account only the first term of the series. It is the the leading term. Let us rewrite l|135p as 



GM 



1 + aie 



-r/Li 



(173) 



The effect of non-Newtonian term can be parameterized by {ai, Li} which could be a useful parameterisation which 
respect to our previous {ai, 02} or {Ge//, L} with Geff = 3G/(4ai). For large distances, where r ^ Li, the 
exponential term vanishes and the gravitational coupHng is G. If r <C I/i, the exponential becomes 1 and, by 
differentiating Eq. (|173p and comparing with the gravitational force measured in laboratory, we get 



Glab — G 



1 + ai 1 



G(l + ai), 



(174) 



where Giab = 6.67 x 10~^ g~^cm'^s~^ is the usual Newton constant measured by Cavendish-like experiments. Of 
course, G and G/a6 coincide in the standard Newtonian gravity. It is worth noticing that, asymptotically, the inverse 
square law holds but the measured coupling constant differs by a factor (1 -|-ai). In general, any correction introduces 
a characteristic length that acts at a certain scale for the self-gravitating systems as in the case of galaxy cluster 
which we are examining here. The range of of the /cth-component of non-Newtonian force can be identified with 
the mass of a pseudo-particle whose effective Compton's length can be defined as 



TOfcC 



(175) 



The interpretation of this fact is that, in the weak energy Hmit, fundamental theories which attempt to unify gravity 
with the other forces introduce, in addition to the massless graviton, particles with mass which also carry the gravi- 
tational interaction [sOl- See, in particular, ^l| for /(i?)-gravity. These masses are related to effective length scales 
which can be parameterized as 



Lk = 2x IQ- 
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There have been several attempts to experimentally constrain Lk and ak (and then m^) by experiments on scales in 
the range 1 cm < r < 1000 km, using different techniques [69, 73, 14^]. In this case, the expected masses of particles 
which should carry the additional gravitational force are in the range 10~^^eV < nik < 10~^ eV. The general outcome 
of these experiments, even retaining only the term k = 1, is that geophysical window between the laboratory and the 
astronomical scales has to be taken into account. In fact, the range 



\ai\ 



Li - 10^^ lO'^ m. 



(177) 



is not excluded at all in this window. An interesting suggestion has been given by Fujii [77||, which proposed that 
the exponential deviation from the Newtonian standard potential could arise from the microscopic interaction which 
couples the nuclear isospin and the baryon number. 

The astrophysical counterparts of these non-Newtonian corrections seemed ruled out till some years ago due to the 
fact that experimental tests of General Relativity seemed to predict the Newtonian potential in the weak energy limit, 
"inside" the Solar System. However, as it has been shown, several alternative theories seem to evade the Solar System 
constraints (see |[31|, and the references therein for recent results) and, furthermore, indications of an anomalous, 
long-range acceleration revealed from the data analysis of Pioneer 10/11, Galileo, and Ulysses spacecrafts (which are 
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Figure 6: Density vs ai: predictions on the behavior of ai. The horizontal black bold line indicates the Newtonian-limit, ai — > 3/4 which 
we expect to be realized on scales comparable with Solar System. Vertical lines indicate typical approximated values of matter density 
(without dark matter) for different gravitational structures: Universe (large dashed) with critical density pcrit ~ 10~^^ kg/m^; galaxy 
clusters (short dashed) with Pd ~ 10"^'^ kg/m'^; galaxies (dot-dashed) with pgai ~ 10~^^ kg/m^; sun (dotted) with psun ~ 10^ kg/m'^. 
Arrows and boxes show the predicted trend for ai. 



now almost outside the Solar System) makes these Yukawa-Hke corrections come again into play [Gj. Besides, it is 
possible to reproduce phenomenologically the flat rotation curves of spiral galaxies considering the values 

ai = -0.92, Li-40kpc. (178) 

The main hypothesis of this approach is that the additional gravitational interaction is carried by some ultra-soft 
boson whose range of mass is nii ~ 10"^'' 10~^^eV. The action of this boson beco mes efficient at galactic scales 
without the request of enormous amounts of dark matter to stabilize the systems [l38l |. 

Furthermore, it is possible to use a combination of two exponential correction terms and give a detailed explanation 
of the kinematics of galaxies and galaxy clusters, again without dark matter model [69j . 

It is worthwhile to note that both the spacecrafts measurements and galactic rotation curves indications come from 
"outside" the usual Solar System boundaries used up to now to test General Relativity. However, the above results do 
not come from any fundamental theory to explain the outcome of Yukawa corrections. In their contexts, these terms 
are phenomenological. 

Another important remark in this direction deserves the fact that some authors [l04l | inte rpre t also the experiments 
on cosmic microwave background like the experiment BOOMERANG and WMAP (HsI . Il49l | in the framework of 
modified Newtonian dynamics again without invoking any dark matter model. 

All these facts point towards the line of thinking that also corrections to the standard gravity have to be seriously 
taken into account beside dark matter searches. 

In our case, the parameters oi 2, which determine the gravitational correction and the gravitational coupling, come 
out "directly" from a field theory with the only requirement that the effective action of gravity could be more general 
than the Hilbert-Einstein theory f{R) — R. This main hypothesis comes from fundamental physics motivations due 
to the fact that any unification scheme or quantum field theory on curved space have to take into account higher 
order terms in curvature invariants [TB|. Besides, several recent results point out that such corrections have a main 
role also at astrophysical and cosmological scales. For a detailed discussion, see [33. It^ . fll^ . 

With this philosophy in mind, we have plotted the trend of oi as a function of the density in FigO As one can see, 
its values are strongly constrained in a narrow region of the parameter space, so that oi can be considered a "tracer" 
for the size of gravitational structures. The value of ai range between {0.8 0.12} for larger clusters and {0.4 0.6} 
for poorer structures (i.e. galaxy groups like MKW4 and RXJ1159). We expect a particular trend when applying 
the model to different gravitational structures. In Fig. [6l we give characteristic values of density which range from 
the biggest structure, the observed Universe (large dashed vertical fine), to the smallest one, the Sun (vertical dotted 
line), through intermediate steps like clusters (vertical short dashed line) and galaxies (vertical dot-dashed line). The 
bold black horizontal line represents the Newtonian limit ai = 3/4 and the boxes indicate the possible values of oi 
that we obtain by applying our theoretical model to different structures. 
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Figure 7: Single temperature fit to tiie total cluster spectrum (upper panel) and total cluster mass within rsoo (given as a fu nctio n of 
M0) (lower panel) are plotted as a function of the characteristic gravitational length L. Temperature and mass values are from [l65l |. 



Similar considerations hold also for the characteristic gravitational length L directly related to both oi and 02. The 
parameter 02 shows a very large range of variation {—10^ — 10} with respect to the density (and the mass) of the 
clusters. The value of L changes with the sizes of gravitational structure (see Fig.[7|, so it can be considered, beside 
the Schwarzschild radius, a sort of additional gravitational radius. Particular care must be taken when considering 
Abell 2390, which shows large cavities in the X-ray surface brightness distribution, and whose central region, highly 
asymmetric, is not expected to be in hydrostat ic eq uilibrium. All results at small and medium radii for this cluster 
could hence be strongly biased by these effects [l65|; the same will hold for the resulting exceptionally high value of 
L. Fig. [7| shows how observational properties of the cluster, which well characterize its gravitational potential (such 
as the average temperature and the total cluster mass within rsoo, plotted in the left and right panel, respectively), 
well correlate with the characteristic gravitational length L. 

For clusters, we can define a gas-density-weighted and a gas-mass-weighted mean, both depending on the series 
parameters ai,2- We have: 



<L>p = 318 kpc 
<L>M = 2738 kpc 



< 02 >p= -3.40 • 10^ 
< 02 >M= -4.15 • 10^ 



(179) 



It is straightforward to note the correlation with the sizes of the cluster cD-dominated-central region and the "gravi- 
tational" interaction length of the whole cluster. In other words, the parameters ai_2, directly related to the first and 
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second derivative of a given analytic /(i?)-model determine the characteristic sizes of the self gravitating structures. 

XIII. WHAT WE HAVE LEARNT FROM CLUSTERS 

We have investigated the possibility that the high observational mass-to-light ratio of galaxy clusters could be 
addressed by f{R)- gravity without assuming huge amounts of dark matter. We point out that this proposal comes 
out from the fact that, up to now, no definitive candidate for dark-matter has been observed at fundamental level 
and then alternative solutions to the problem should be viable. Furthermore, several results in /(i?)-gravity seem to 
confirm that valid alternatives to ACDM can be achieved in cosmology. Besides, as discussed in the Introduction, the 
rotation curves of spiral galaxies can be explained in the weak field limit of /(i?)-gravity. Results of our analysis go 
in this direction. 

We have chosen a sample of relaxed galaxy clusters for which accurate spectroscopic temperature measurements and 
gas mass profiles are available. For the sake of simplicity, and considered the sample at our disposal, every cluster has 
been modelled as a self-bound gravitational system with spherical symmetry and in hydrostatic equihbrium. The mass 
distribution has been described by a corrected gravitational potential obtained from a generic analytic /(i?)-theory. 
In fact, as soon as f{R) ^ R, Yukawa-like exponential corrections emerge in the weak field limit while the standard 
Newtonian potential is recovered only for f{R) = R, the Hilbert-Einstein theory. 

Our goal has been to analyze if the dark-matter content of clusters can be addressed by these correction potential 
terms. As discussed in detail in the previous sections and how it is possible to see by a rapid inspection of figures, the 
clusters of the sample are consistent with the proposed model at Icr confidence level. This shows, at least qualitatively, 
that the high mass-to-Hght ratio of clusters can be explained by using a modified gravitational potential. The good 
agreement is achieved on distance scales starting from 150 kpc up to 1000 kpc. The differences observed at smaller 
scales can be ascribed to non-gravitational phenomena, such as cooling fiows, or to the fact that the gas mass is not 
a good tracer at this scales. The remarkable result is that we have obtained a consistent agreement with data only 
using the corrected gravitational potential in a large range of radii. In order to put in evidence this trend, we have 
plotted the baryonic mass vs radii considering, for each cluster, the scale where the trend is clearly evident. 

In our knowledge, the fact tha t .f (i?)-g ravity could work at these scales has been only supposed but never achieved 
by a direct fitting with data (see [lp,llO0| for a review). Starting from the series coefficients oi and 02, it is possible to 
state that, at cluster scales, two characteristic sizes emerge from the weak field limit of the theory. However, at smaller 
scales, e.g. Solar System scales, standard Newtonian gravity has to be dominant in agreement with observations and 
experiments. 

In summary, if our considerations are right, gravitational interaction depends on the scale and the infrared limit is 
led by the series coefficient of the considered effective gravitational Lagrangian. Roughly speaking, we expect that 
starting from cluster scale to galaxy scale, and then down to smaller scales as Solar System or Earth, the terms of 
the series lead the clustering of self-gravitating systems beside other non-gravitational phenomena. In our case, the 
Newtonian limit is recovered for ai — > 3/4 and L{ai,a2) » r at small scales and for L{ai,a2) ^ r at large scales. 
In the first case, the gravitational coupHng has to be redefined, in the second Goo — G. In th ese l imits, the Hnear 
Ricci term is dominant in the gravitational Lagrangian and the Newtonian gravity is restored |l29l |. Reversing the 
argument, this could be the starting point to achieve a theory capable of explaining the strong segregation in masses 
and sizes of gravitationally-bound systems. 

XIV. CONCLUSIONS 

The present status of art of cosmology shows that the Standard Cosmological Model, based on General Relativity, 
nucleosynthesis, cosmic abundances and large scale structure, has some evident difficulties. These ones, first of all, 
rely on some lack of a self-consistent formulation of missing matter and cosmic acceleration issues; such shortcomings 
give rise to further difficulties in interpreting observational data. With an aphorism, one can say that we have a book, 
but not the alphabet to read it. 

Nowadays there two main philosophical approaches aimed to solve this problem. From one side, there are researchers 
which try to solve shortcomings of Standard Cosmological Model assuming that General Relativity is right but we 
need some exotic, invisible kinds of energy and matter to explain cosmic dynamics and large scale structure. On the 
other side, there are people which beheve that General Relativity is not the definitive and comprehensive theory of 
gravity, and that it should be revised at ultraviolet scales (quantum gravity) and infrared scales (extragalactic and 
cosmic scales). In the latter case, dark energy and dark matter could be nothing else but the signals that we need 
a more general theory at large scales, also if General Relativity works very well up to Solar System scales. To some 



39 



5x10 



1x10 
--5 X 10 



Si X 10 
5x10 



1x10 




100 150 200 300 500 700 1000 
r (kpc) 



Figure 8: As an example of the above results, we have plotted the baryonic mass vs radii for Abell A133. Dashed line is the experimental- 
observed estimation Eq. H164|l of baryonic matter component (i.e. gas, galaxies and cD-galaxy); solid line is the theoretical estimation 
Eq. II163II for baryonic matter component. Dotted lines are the l-cr confidence levels given by errors on fitting parameters plus statistical 
errors on mass profiles as discussed in $ IXI Dl in the right panel. 
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Figure 9: As the above case, for cluster Abell 383. 



extent, this could be seen as a sort of philosophical debate without solution, but there are possibilities to move the 
question toward a physical viewpoint. 

The /(i?)-gravity is strictly related to the second point of view. It is a fruitful approach to generalize General 
Relativity towards the solution also if, most of the models in literature are nothing else but phenomenological models. 
It is interesting to note that as soon as Einstein formulated his General Relativity, many authors (and Einstein himself) 
started to explore other possibilities (see [33] for a review). At the beginning, these researches were mainly devoted to 
check the mathematical consistency of General Relativity but the issues to achieve the unification of gravity with the 
other interactions (e.g. electromagnetism) pushed several authors to develop alternative gravity theories. Today, one 
of the goals of alternative gravity is to understand the effective content and dynamics of the Universe. This question 
is recently become dramatic since assuming that more than 95% of cosmic matter-energy is unknown at fundamental 
level is highly disturbing. Alternative gravity could be a way out to this situation. The present status of observations, 
also if we are living in the era of Precision Cosmology, does not allow in discriminating between alternative gravity, 
from one side, and the presence of dark energy and dark matter, from the other side (the forthcoming LHC experiments 
should aid in this sense if new fundamental particles will be detected) . 

However, as discussed in this review, cosmography may be a useful tool to discriminate among different cosmological 
models being, by definition, a model-independent approach: any cosmographic parameter can be estimated without 
assigning an a priori cosmological model. So cosmography can be used in two ways: 

• One can use it to discriminate between General Relativity and alternative theories. This issue strictly depends on 
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Figure 10: As the above cases, for cluster Abell 478. 




the possibility to have good quahty data at disposal. We need some minimum sensibility and error requirements 
on data surveys to solve this question. At the moment, we have not them and we are not able to do this since 
standard candles are not available at very high red shifts [33 |. 

• We can use the cosmographic parameters to constraint cosmological models as we have done in this paper for 
/(i?)-gravity. Being these parameters model-independent, they results natural "priors" to any theory. As above, 
the accuracy in estimating them is a crucial issue. 

We have used, essentially, SNela but other classes of objects have to be considered in order to improve such an 
accuracy (e.g. CMBR, bright galaxies, GRBs, BAOs, weak lensing and so on). Forthcoming space missions will be 
extremely useful in this sense. 

Beside cosmography, we have discussed also if /(i?)-gravity could be useful to address the problem of mass profile 
and dynamics of galaxy clusters. This issue is crucial in view of achieving any correct model for large scale structure. 

Taking into account the weak-field limit of a generic analytic /(i?)-function, it is possible to obtain a scale- dependent 
gravity, where scales of self-gravitating systems could naturally emerge. In this way, one could successfully explain 
dark matter profiles ranging from galaxies to clusters of galaxies. The results are preHminary but seems to indicate 
a way in which the dark matter puzzle could be completely solved. 

In conclusion, the main lesson of this work is that since it is very difficult to discriminate among the huge amount of 
cosmological models which try to explain the data (deductive approach), it could be greatly fruitful to "reconstruct" 
the final cosmological model by an inductive approach, that is without imposing it a priori but adopting the philosophy 
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Figure 12: As the above cases, for cluster Abell 2029. 



to use the minimum number of parameters This "inverse scattering approach" could be not fully satisfactory but 
could lead to self-consistent results. 
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